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1?  ABSTRACT  20"  .sees: 

This  report  documents  the  development  and  testing  of  an  advanced  geophysical 
parameter  retrieval,  analysis  and  prediction  system  for  global  numerical  weather 
prediction.  The  additions  and  improvements  to  the  existing  data  assimilation 
system  at  the  Phillips  Laboratory  are  the  following:  a  new  surface  pressure 
analysis,  modifications  for  high  resolution  analysis,  an  optimal  retrieval  of 
SSM/T  data  based  on  Bayesian  estimation,  use  of  precipitation  data  to  improve 
the  ireaisture  and  divergence  analysis  and  as  input  to  a  diabatic  nonlinear  normal 
mode  initialization.  The  system  uses  a  version  of  the  advanced  physics  global 
spectral  model  of  PL. 

The  results  of  simulation  tests  show  a  significant  improvement  of  the  moisture 
analyses,  especially  in  the  tropics  and  Southern  Hemisphere,  and  a  smaller,  but 
detectable  improvement  of  the  other  variables,  except  in  the  stratosphere. 
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Fig.  1:  Data  flow  diagram  for  CASCADE.  The  black  arrows  represent  the 
model  variables,  the  white  arrows  the  data,  and  the  gray  arrows 
auxiliary  variables  that  are  computed  by  the  system 

Fig.  2:  Impact  of  improved  surface  pressure  analysis  on  500  hPa  rms  analysis 
error  for  the  Northern  Hemisphere  extratropics.  The  surface  pressure 
experiment  (dashed  line)  starts  on  Julian  day  325,  from  a  state  identical 
to  the  control  experiment  (solid).  The  RAOB  observation  error 
standard  deviation  (dotted)  is  shown  for  reference. 

Fig.  3:  Example  of  retrieval  for  a  clear  case.  The  solid  lines  are  the  first  guess 
errors  of  temperature  (left)  and  logarithm  of  specific  humidity  (right). 
The  dashed  lines  are  the  corresponding  retrieval  errors 

Fig.  4:  Same  as  Fig.  3,  but  for  a  cloudy  case. 

Fig.  5:  Schematic  showing  the  sequence  of  spinup  forecast,  spinup 

assimilation,  7-day  assimilation,  and  4-day  forecasts  from  day  3,  5,  and 
7. 

Fig.  6:  Root  mean  square  errors  of  500  hPa  geopotential  height  computed  for 
the  globe  (a).  Northern  Hemisphere  extra  tropics  (b).  Tropics  (c),  and 
Southern  Hemisphere  extratropics  (d).  The  time  period  covers  the 
spinup  forecast,  spinup  assimilation,  the  7-day  assimilation  and  the  4- 
day  forecasts  generated  from  it. 

Fig.  7:  Root  mean  square  errors  of  700  hPa  relative  humidity  computed  for  the 
globe  (a).  Northern  Hemisphere  extra  tropics  (b),  Tropics  (c),  and 
Southern  Hemisphere  extratropics  (d).  The  time  period  covers  the 
spinup  forecast,  spinup  assimilation,  the  7-day  assimilation  and  the  4- 
day  forecasts  generated  from  it. 

Fig.  8:  Root  mean  square  errors  of  700  hPa  horizontal  wind  computed  for  the 
globe  (a).  Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and 
Southern  Hemisphere  extratropics  (d).  The  time  period  covers  the 
spinup  forecast,  spinup  assimilation,  the  7-day  assimilation  and  the  4- 
day  forecasts  generated  from  it. 

Fig.  9:  Root  mean  square  errors  of  300  hPa  horizontal  wind  computed  for  the 
globe  (a).  Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and 
Southern  Hemisphere  extra  tropics  (d).  The  time  period  covers  the 
spinup  forecast,  spinup  assimilation,  the  7-day  assimilation  and  the  4- 
day  forecasts  generated  from  it. 
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Fig.  10:  Vertical  profiles  of  the  root  mean  square  errors  of  geopotential  height 
averaged  over  the  last  5  days  of  the  assimilation,  for  the  globe  (a). 
Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and  Southern 
Hemisphere  extratropics  (d). 

Fig.  11:  Vertical  profiles  of  the  root  mean  square  errors  of  temperature 

averaged  over  the  last  5  days  of  the  assimilation,  for  the  globe  (a). 
Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and  Southern 
Hemisphere  extra  tropics  (d). 

Fig.  12:  Vertical  profiles  of  the  root  mean  square  errors  of  relative  humidity 
averaged  over  the  last  5  days  of  the  assimilation,  for  the  globe  (a). 
Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and  Southern 
Hemisphere  extra  tropics  (d). 

Fig.  13:  Vertical  profiles  of  the  root  mean  square  errors  of  the  horizontal  wind 
(vector  difference)  averaged  over  the  last  5  days  of  the  assimilation,  for 
the  globe  (a).  Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and 
Southern  Hemisphere  extratropics  (d). 

Fig.  14:  Root  mean  square  errors  of  500  hPa  geopotential  height  as  a  function 
of  forecast  lead  time,  averaged  over  the  day  3, 5,  and  7  forecasts,  for  the 
globe  (a).  Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and 
Southern  Hemisphere  extratropics  (d). 

Fig.  15:  Root  mean  square  errors  of  850  hPa  relative  humidity  as  a  function  of 
forecast  lead  time,  averaged  over  the  day  3,  5,  and  7  forecasts,  for  the 
globe  (a).  Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and 
Southern  Hemisphere  extratropics  (d). 

Fig.  16:  Root  mean  square  errors  of  850  hPa  zonal  wind  as  a  function  of 

forecast  lead  time,  averaged  over  the  day  3,  5,  and  7  forecasts,  for  the 
globe  (a).  Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and 
Southern  Hemisphere  extratropics  (d). 

Fig.  17:  Root  mean  square  errors  of  300  hPa  zonal  wind  as  a  function  of 

forecast  lead  time,  averaged  over  the  day  3,  5,  and  7  forecasts,  for  the 
globe  (a).  Northern  Hemisphere  extratropics  (b),  Tropics  (c),  and 
Southern  Hemisphere  extratropics  (d). 

Fig.  A.l:  Distribution  of  CDW  pressures. 

Fig.  A. 2:  Horizontal  correlations  for  low  level  (1000-600  hPa)  CDWs.  Shown 
are  correlations  for  the  background  field  evaluated  at  the  CDW 
locations  (top),  original  CDW  errors  (middle)  and  filtered  CDW  errors 
(bottom)  for  u  (left)  and  v  (right)  wind  components. 
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Fig.  A.3:  Horizontal  correlations  for  high  level  (600-100  hPa)  CDWs.  Same  as 
Fig.  A.2. 

Fig.  B.l:  Mean  (a,  c)  and  rms  (b,  d)  errors  of  temperature  (a,  b)  and  relative 
humidity  (c,  d)  before  and  after  the  retrieval  (dotted  and  solid  lines, 
respectively).  Statistics  are  from  the  7-day  assimilation  for  all  analysis 
grid  points  for  which  retrievals  were  computed. 

Fig  C.l:  Profiles  of  divergence  along  o  surfaces  for  5  different  categories  of 

rainfall  rate  for  grid  points  in  the  tropics.  Category  1  is  shown  in  solid, 
2  in  dotted,  3  in  dashed,  4  in  dash-dotted,  and  5  in  dash-dot-dot-dot. 


Table  C.l:  Definition  of  categories  1-5.  Riow  and  Rupp  are  the  lower  and  upper 
limits  of  rain  rate  (note  that  10'®  m/s  corresponds  to  0.864  mm/day). 

Table  C.2:  Configuration  of  initialization  experiments.  Convective  heating  and 
precipitation  scaled  by  ratio  of  observed  to  predicted  precipitation  during 
NMI  and  forecast  in  all  experiments  other  than  CONTROL. 

Temperature  and  moisture  adjustment  refers  to  the  adjustment  of 
lowest  layer  values  only. 


1  Introduction 


This  report  documents  the  development  and  testing  of  an  advanced 
geophysical  parameter  retrieval,  analysis  and  prediction  system  for  global 
numerical  weather  prediction  (NWP),  performed  under  the  Phillips 
Laboratory  contract  No.  F19628-89-C-0044. 

We  called  this  system  "CASCADE",  because  of  the  shape  of  the  data  flow 
diagram  (see  Figure  1)  and  because  the  work  was  designed  to  provide  a  series 
of  incremental  improvements  to  the  system,  each  flowing  into  the  next  one. 
CASCADE  is  a  step  towards  an  advanced  and  nearly  ideal  analysis  system  for 
global  NWP  which,  we  feel,  is  attainable  by  the  mid  to  late  1990's. 

The  modular  design  of  CASCADE  provided  a  basis  for  our  work,  which 
allowed  us  to  make  use  of  previous  efforts  and  to  concentrate  on  the  aspects 
of  the  analysis  system  that  are  most  important  to  the  Air  Force  concerns  and 
that  may  not  be  fully  addressed  by  the  op)erational  meteorological  centers. 
Modularity  implies  that  the  overall  procedure  can  be  easily  reconfigured  by 
replacing  sub-procedures.  CASCADE  was  developed  step  by  step, 
incorporating: 

•  a  variable  resolution  optimal  interpolation  scheme; 

•  an  enhanced  satellite  data  handling  capability  based  on  an  optimal 
retrieval  approach  to  process  the  DMSP  sensor  suite  data  including  the 
SSM/T,  SSM/T-2,  and  eventually  SSM/I  and  OLS; 

•  algorithms  to  use  satellite  observed  clouds  to  infer  precipitation  rates. 
The  inferred  precipitation  rates  are  used  to  deduce  latent  heating  profiles  and 
to  adjust  the  relative  humidity  field  in  convective  regions; 

•  inclusion  and  testing  of  diabatic  effects  in  the  nonlinear  normal  mode 
initialization. 

Four  simulated  data  assimilation  experiments  were  performed  to  test  the 
various  elements  of  the  CASCADE  system.  These  experiments  were 
performed  with  data  extracted  from  a  "Nature  run"  provided  by  the  European 
Centre  for  Medium-Range  Weather  Forecasts  (ECMWF). 
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2 _ Background 

Remote  sensing  from  space  and  numerical  weather  prediction  technology 
and  practice  have  evolved  significantly  since  their  operational  beginnings  in 
the  1950s  and  1960s.  Inferences  of  atmospheric  temperature  and  moisture, 
clouds,  precipitation,  winds,  and  surface  properties  can  now  be  made  from  a 
variety  of  visible,  infrared,  and  microwave  satellite  sensors.  For  global  NWP 
it  is  primarily  the  retrieved  temperature  profile  and  wind  data  which  have 
been  used.  The  usefulness  of  other  geophysical  parameters  for  NWP  is  not  yet 
well  established  and  is  largely  untested,  although  a  great  many  retrieval 
methods  have  been  proposed  or  developed  for  a  variety  of  potentially 
interesting  parameters.  The  reason  for  this  is  twofold.  On  the  one  hand,  these 
data  are  expected  to  have  the  most  impact  on  mesoscale  modeling,  but 
mesoscale  analysis  systems  are  still  relatively  primitive.  On  the  other  hand, 
moisture  variables  -  that  is,  specific  humidity,  clouds  and  precipitation  -  are 
retrievable  and  are  potentially  very  useful,  but  are  not  easily  assimilated  by 
current  methods.  It  is  theoretically  possible  to  retrieve  specific  humidity 
profiles  by  using  methods  analogous  to  those  used  to  retrieve  temperatures. 
However,  results  to  date  with  available  sensors  have  not  been  wholly 
satisfactory. 

Other  moisture  parameters  are  somewhat  easier  to  retrieve;  for  example, 
cloud  cover  may  be  obtained  by  a  variety  of  techniques  from  imagery  in 
visible  and/or  infrared  spectral  regions,  and  vertically  integrated  water  vapor 
is  well  correlated  with  emissions  at  microwave  frequencies,  especially  over 
the  ocean.  Unfortunately,  these  parameters  are  not  readily  assimilated  by  a 
global  NWP  model,  and  special  methods  are  required.  A  variety  of  other 
meteorological  parameters  which  may  be  used  for  NWP  are  retrievable. 

These  include  winds  at  the  surface  from  microwave  sensors  and  cloud  drift 
winds  aloft,  as  well  as  other  surface  properties,  such  as  soil  moisture,  albedo, 
snow  cover,  temperature,  and  fluxes.  The  impact  of  current  observing 
systems  on  global  NWP  is  not  large  except  in  the  southern  hemisphere, 
where  conventional  data  sources  are  limited.  Increasing  the  accuracy  and 
vertical  resolution  of  the  sensors  and  improving  the  retrieval  and 
assimilation  systems  should  lead  to  better  forecasts.  We  believe  that  only  a 
small  part  of  the  potential  usefulness  of  satellite  remote  sensing  for  global 
NWP  has  been  realized. 

It  is  evident  that  improvements  in  NWP  come  from  parallel  progress  in  the 
forecast  models,  their  resolution,  numerical  schemes  and  simulation  of  the 
physics,  and  in  the  analysis  of  the  data  that  provide  the  initial  conditions  for 
the  forecasts.  Data  analysis  techniques  have  evolved  considerably  over  the 
last  few  years,  partly  because  of  the  introduction  of  new  data  types  such  as 
satellite  data.  The  development  of  the  optimal  interpolation  (OI)  method,  for 
example,  was  spurred  by  the  fact  that  different  data  typos  have  different  error 
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characteristics.  Nevertheless  most,  if  not  all,  of  the  existing  operational  data 
analysis  methods  are  still  closely  tied  to  the  conventional  type  of  data,  namely 
radiosoundings,  which  provide  profiles  of  the  variables  used  by  the  NWP 
models,  simultaneously  at  many  points  on  the  globe.  Consequently  much 
effort  has  been  devoted  in  trying  to  make  all  data  look  like  radiosoundings: 
the  majority  of  all  satellite  data  retrieval  studies  have  been  devoted  to  this 
goal. 

It  is  being  realized  more  and  more,  however,  that  trying  to  make  all  data  look 
like  radiosoundings  may  not  be  the  best  tactic.  Indeed  new  types  of  data  have 
fundamentally  different  characteristics.  The  most  impc  tant  differences  are 
that  they  are  asynoptic,  nonlocal,  indirect,  and  often  incomplete  (e.g.  Doppler 
radar).  Instead  of  being  taken  all  at  the  same  time  like  the  synoptic  data,  they 
are  taken  more  or  less  continuously  or  at  random  times.  Remote  sensing 
techniques  sample  a  much  larger  volume  of  the  atmosphere  than  do 
radiosondes,  and  the  radiances  that  are  measured  are  only  indirectly  related  to 
the  variables  needed  by  the  NWP  models,  through  a  convolution  of  the 
profiles  by  a  complicated  transfer  function.  In  some  cases  (e.g.  imagery  data), 
very  little  NWP  use  has  been  made  of  vast  quantities  of  data  because  of  the 
difficulty  of  making  them  look  like  a  radiosounding.  The  current  thinking  is 
that  it  would  be  preferable  to  use  data  directly,  in  the  form  and  at  the  time 
they  are  taken.  Different  approaches  have  been  taken,  mainly  using 
variational  or  Bayesian  methods. 

Another  important  problem  of  NWP  is  the  analysis  of  moisture  data.  The 
humidity  field  tends  to  have  a  high  variability  at  small  scales,  so  that  its 
determination  from  conventional  data  is  rather  poor.  Up  to  now,  not  much 
use  has  been  made  of  satellite  data  to  help  specify  the  moisture  field.  The 
result  is  a  poor  humidity  analysis  and,  consequently,  a  poor  humidity 
forecast.  The  first  few  hours  of  a  forecast  are  usually  characterized  by  a  rapid 
re-adjustment  of  the  moisture  field  to  one  compatible  with  the  model 
thermodynamic  structure.  This  process  is  accompanied  by  large  variations  in 
the  precipitation  rates  and  large  changes  in  the  diagnosed  cloud  cover.  This  is 
part  of  the  so-called  "spin-up"  problem. 

Related  to  this  problem  is  the  question  of  properly  initializing  or  balancing 
the  model  mass  and  wind  fields  to  avoid  the  development  of  large  amplitude 
gravity  waves.  In  this  balance,  one  needs  to  take  into  account  the  effect  of 
diabatic  heating,  especially  the  latent  heat  release.  It  is  obvious  that,  with  a 
poor  moisture  analysis,  the  diabatic  heating  cannot  be  computed  properly. 
Conversely,  if  t.e  latent  heat  release  can  be  estimated  (e.g.  from  imagery  data), 
then  a  better  moisture  analysis  may  be  made,  more  in  balance  with  the 
thermodynamic  field,  thus  reducing  the  spin-up  problem. 

The  ultimate  goal  of  data  retrieval  and  analysis  is  to  find  a  way  of  combining 
all  the  information  that  we  have  about  the  atmosphere  and  its  processes  in  a 


3 


consistent  manner  and  properly  accounting  for  all  aspects  of  the  data.  The 
methodology  to  achieve  this  goal  will  most  likely  be  a  four-dimensional  (4-D) 
variational  analysis  of  all  available  data,  using  an  adjoint  model  to  determine 
the  gradient  of  the  objective  function  (measuring  the  closeness  between  the 
model  state  and  the  measurements)  with  respect  to  the  initial  state.  The  data 
would  include  imagery,  radiances  and  all  conventional  data.  New  data 
sources  such  as  ground  based  lidar  profilers  with  fine  temporal  resolution, 
space  borne  Doppler  lidar  wind  profilers,  etc.  would  be  accommodated  with 
ease.  Within  such  a  scheme,  improvements  to  the  model  immediately  yields 
improvements  to  the  retrieval /analysis. 

Most  operational  centers,  including  the  National  Meteorological  Center 
(NMC),  ECMWF  and  the  French  weather  service,  are  working  on  this  in  a 
straightforward  manner,  developing  the  adjoint  model  methodology 
mentioned  above.  Although  developing  such  a  methodology  in  a  few  years 
was  an  effort  that  by  far  exceeded  the  resources  available  to  us,  we  have  kept 
this  ultimate  analysis  in  mind  in  our  research  effort.  The  work  presented 
here  on  the  development  of  CASCADE  represents  several  steps  toward  this 
ultimate  objective,  making  the  best  possible  use  of  the  work  already 
completed  by  AER  and  Phillips  Laboratory  personnel,  and  keeping  in  focus 
the  specific  concerns  of  the  US  Air  Force. 

In  the  design  of  CASCADE  we  have  followed  some  general  principles  which, 
we  feel,  are  essential  for  the  development  of  the  ideal  data  assimilation 
method  of  the  future.  The  main  points  are  the  following: 

(1)  Use  of  all  available  data.  Although  using  all  data  is  obviously  desirable,  we 
have  seen  that,  for  various  reasons,  many  types  of  data  are  currently  not  used 
at  all.  CASCADE  does  not  achieve  this  requirement  either,  but  it  uses  more 
fully  the  data  from  the  DMSP  suite  of  instruments,  as  well  as  imagery  data  (or 
RTNEPH  data  as  a  surrogate  for  imagery). 

(2)  Use  of  remotely  sensed  physical  forcing  mechanisms  such  as  latent 
heating.  For  example,  since  the  effect  of  latent  heat  release  is  important  in 
intense  storms,  it  is  natural  to  consider  using  observed  precipitation  rates  to 
force  the  model  by  an  externally  specified  heating  function  during  an 
assimilation  period.  In  CASCADE,  a  first  step  is  taken  by  using  estimates  of 
latent  heating  that  can  be  derived  from  remotely  sensed  data  in  the  nonlinear 
normal  mode  initialization. 

(3)  Direct  use  of  remotely  sensed  data  through  a  variational  analysis  rather 
than  by  a  retrieval  method.  Such  a  procedure  avoids  the  introduction  of 
errors  by  the  retrieval  technique  and  allows  optimal  use  (in  a  Bayesian  sense) 
of  the  data.  This  has  been  the  main  thrust  of  our  effort. 
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While  CASCADE  is  a  step  towards  an  ideal  data  assimilation  system,  it  still 
falls  short  in  a  number  of  ways.  The  most  important  principles  that 
CASCADE  does  not  follow  are: 

(1)  Using  data  continuously  without  time  lags.  The  conventional  practice  of 
batching  data  into  6  h  increments  is  still  followed  by  CASCADE.  This  is 
reasonable  as  long  as  only  the  largest  scales  are  under  consideration.  With 
increasingly  fine  mesh  models  and  high  resolution  observations,  smaller 
scales  will  become  more  important,  and  a  6  hour  data  increment  will  blur  the 
observational  picture. 

(2)  Reducing  the  need  for  normal  mode  initialization.  Normal  mode 
initialization  procedures  are  reasonable  approaches  when  certain 
assumptions  are  met,  but  these  assumptions  are  less  valid  in  the  tropics 
where  the  Rossby  number  is  not  small,  and  in  the  boundary  layer  where  the 
effect  of  friction  is  strong.  Nevertheless,  the  design  of  CASCADE  still  requires 
the  use  of  an  initialization. 

These  two  unmet  requirements  can  be  resolved  only  with  a  fully 
quadridimensional  data  assimilation  system.  For  the  time  being,  CASCADE 
remains  a  tridimensional  scheme,  doing  data  analysis  intermittently  in  six 
hour  intervals. 
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3  Experimental  Design 

3.1  The  design  of  CASCADE 

CASCADE  is  described  in  block  form  by  Fig.  1.  This  is  a  data  flow  description. 
The  arrows  in  the  figure  represent  data  being  transferred;  black  when  the  data 
are  model  states  (i.e.  a  complete  description  of  the  model  atmosphere: 
normally,  u,  v,  Z,  p*,  rh),  white  when  the  data  are  observations  and  gray 
when  the  data  are  auxiliary  variables  that  are  derived  either  from  the  model 
variables  or  the  observations.  Each  step  of  CASCADE  adds  new  (i.e. 
independent)  information.  Wherever  possible,  the  background  field  and 
observations  are  combined  to  minimize  the  probable  analysis  error. 


Figure  1:  Data  flow  diagram  for  CASCADE.  The  black  arrows  represent 
the  model  variables,  the  white  arrows  the  data,  and  the  gray  arrows 
auxiliary  variables  that  are  computed  by  the  system. 

CASCADE  starts  with  a  standard  optimal  interpolation  analysis  of  the 
conventional  data  and  other  usual  data  including  aircraft  data,  cloud  drift 
winds  and  operational  NOAA  TOYS  satellite  temperature  soundings.  At  the 
next  step  (radiance  retrieval),  data  from  the  SSM/T-1  and  SSM/T-2  are  used  to 
improve  the  analysis  of  temperature  and  moisture.  Changes  in  the 
extratropical  temperature  fields  are  reflected  in  the  wind  analysis  by 
geostrophy.  Horizontal  consistency  is  insured  by  an  appropriate  filter.  After 
this  step  cloud  information  is  used  to  refine  the  moisture  analysis.  It  is  also 
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used  to  determine  latent  heating  rates  which  become  a  constraint  in  the 
initialization  that  follows. 

CASCADE  differs  from  an  ideal  scheme  in  two  important  ways.  First  of  all,  as 
has  already  been  mentioned,  it  is  an  intermittent  scheme,  with  all  data 
grouped  in  batches  and  analyzed  together  every  6  hours.  Secondly,  each 
retrieval  is  split  into  a  vertical  (retrieval  proper)  and  horizontal  (smoothing) 
component.  In  an  ideal  system  the  analysis  would  be  a  continuous  process, 
and  the  vertical  and  horizontal  components  of  the  retrievals  would  be 
combined.  Note  that  we  avoid  the  problem  of  retrieval  -  first  guess 
correlations  by  processing  the  different  data  in  series  instead  of  in  parallel. 

CASCADE  was  developed  in  a  series  of  steps,  each  one  adding  a  new  module. 
The  baseline  version  consisted  of  modules  that  were,  for  the  most  part, 
already  available,  although  we  had  to  make  fairly  substantial  improvements 
to  the  existing  OI  to  include  a  surface  pressure  analysis  and  to  use  a  higher 
resolution  model. 

In  the  following  sections  we  present  the  various  modules  making  up 
CASCADE.  In  each  section  we  first  briefly  describe  what  was  already  available 
to  build  up  version  1.0,  then  we  detail  improvements  that  were  made. 

3.1.1  Optimal  interpolation 

The  starting  point  for  CASCADE  was  the  ASAP  optimal  interpolation  scheme 
that  has  been  developed  at  AFGL  (now  the  Phillips  Laboratory).  The  ASAP  OI 
was  developed  by  SASC  for  AFGL/LYP  and  is  documented  in  a  series  of 
reports  by  Norquist  and  others  (Halberstam  et  ai,  1984;  Norquist,  1982; 
Norquist,  1983;  Norquist,  1984;  Norquist,  1986).  The  ASAP  OI  is  based  on  the 
NMC  assimilation  system  (Bergman,  1979;  McPherson  et  al.,  1979),  but  was 
thoroughly  redesigned  and  recoded  by  SASC  personnel.  It  is  a  multivariate 
analysis  of  height  and  wind  components  and  a  univariate  analysis  of  relative 
humidity,  both  in  model  sigma  layers. 

Data  used  by  the  height-wind  analysis  include  Type  1  observations 
(radiosondes,  pibals,  etc.).  Type  2  observations  (aircraft).  Type  4  observations 
(satellite-retrieved  temperatures  or  thicknesses)  and  Type  6  observations 
(cloud  drift  winds  (CDWs)).  The  Type  3  surface  observations  are  not  used  at 
all.  The  moisture  analysis,  in  the  original  ASAP,  uses  only  Type  1  data.  The 
CDW  data  (Type  6  data)  are  combined  (i.e.,  locally  averaged)  into  super-obs. 
There  are  two  principal  reasons  for  doing  this:  first,  to  limit  the  total  number 
of  observations,  so  that  computer  memory  restrictions  are  not  exceeded,  and 
second,  the  CDW  errors  are  strongly  related  horizontally  because  the  main 
error  source  is  due  to  height  assignment  (McPherson,  1984). 

The  version  of  the  OI  used  at  the  start  of  the  project  included 

1)  a  strict  decoupling  of  height  and  wind  analysis  in  the  equatorial  belt. 
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2)  the  simultaneous  solution  for  heights  and  winds  so  increments  are 
geostrophically  balanced  (normal  equations  are  now  solved  using 
LINPAK  routines), 

3)  an  underground  analysis  of  heights  and  winds,  to  assure  smooth 
extrapolation  of  height  analysis  increments  to  the  surface, 

4)  correct  specification  of  the  right  hand  side  of  the  normal  equations,  and 

5)  a  more  efficient  and  accurate  method  of  converting  satellite  thicknesses 
to  level  heights. 

Three  sets  of  modifications  were  made  to  the  OI  code  to  produce  the  baseline 
version:  implementation  of  variable  resolution,  introduction  of  a  surface 
pressure  analysis,  and  modified  data  selection  for  the  lowest  layers. 

3.1. 1.1  Variable  resolution  code  changes 

A  series  of  changes  were  made  to  the  analysis  codes  (ASAP  and  MASAP  for 
moisture)  and  their  associated  pre-  and  postprocessors  to  allow  for  variable 
resolution. 

The  changes  to  the  analysis  codes  were  made  as  a  series  of  small  changes, 
which  enabled  us  to  reproduce  the  results  of  the  initial  code  with  small 
errors,  the  source  of  which  we  clearly  identified.  Significant  changes  have 
been  reflected  in  the  new  documentation  for  ASAP,  which  was  published  as  a 
scientific  report  (Hoffman  et  al,  1988).  The  main  changes  are  the  following; 

All  common  blocks  are  now  incorporated  in  the  code  by  "include" 
statements  and  are  the  same  whenever  they  are  used. 

The  constants  for  LOWTMP  (the  Flattery  algorithm)  are  calculated  once  at 
the  start  of  ASAP,  since  these  constants  depend  on  the  vertical 
structure. 

The  first  guess  error  growth  rates  are  stored  at  mandatory  levels  and 
interpolated  to  the  sigma  structure  instead  of  being  hard-wired  for  a 
particular  model  structure. 

The  code  for  the  3DNeph  and  surface  weather  options  (used  by  Norquist, 
1988)  has  been  deleted. 

The  code  is  now  strictly  standard  FORTRAN. 

Modifications  for  MASAP  paralleled  those  for  ASAP.  The  largest  differences 
during  this  series  of  modifications  are  due  to  the  calculation  of  the  LOWTMP 
constants  to  full  machine  precision  instead  of  specifying  them  on  a  DATA 
statement  to  9  or  10  digits,  and  to  the  change  in  the  error  growth  rate 
computation,  which  produced  differences  detectable  at  the  0(1 g/kg)  level. 

The  pre-  and  postprocessors  (PPPs)  were  completely  redesigned,  in  a  modular 
fashion.  Since  the  global  spectral  model  (GSM)  PPPs  already  exist  in  variable 
resolution  format,  we  only  implemented  our  design  for  the  PPPs  for  the  OI. 
Many  of  the  subroutines  needed  to  affect  the  implementation  for  the  GSM  are 
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now  available,  and  the  two  sets  of  PPPs  could  be  unified  at  some  later  date. 

This  has  been  done  for  the  GSM  postprocessor.  The  subroutines  were  tested 
extensively. 

Three  new  main  modules  were  coded  and  tested: 

•  Spectral  to  grid  p)oint  (OISTG):  Reads  in  spectral  coefficient  data  set  and 
writes  out  the  grid  point  first  guess  files. 

•  Grid  point  to  spectral  (OIGTS):  Reads  the  analysis  grid  point  files,  updates 
the  pressure  field  and  outputs  the  spectral  coefficient  analysis  file. 

•  Spectral  filter  for  estimated  analysis  error  (PREAES).  Reformats  height  and 
relative  humidity  estimated  analysis  errors. 

3.1. 1.2  Surface  pressure  analysis 

Early  in  the  project  we  decided  to  re-examine  the  question  of  how  to  arrive  at 
sigma-layer  temperature  increments  from  an  analysis  of  geopotential  height. 
Results  from  our  previous  contracts  (Hoffman  et  al.,  1990;  Louis  et  ai,  1989) 
had  indicated  that  the  lack  of  a  surface  pressure  analysis,  the  neglect  of  surface 
data,  and  the  consequent  anchoring  of  the  satellite  thicknesses  to  the  first 
guess  surface  pressure  severely  limited  the  applicability  of  the  results  from  our 
OSEs  and  OSSEs.  We  had  been  especially  concerned  about  the  relatively  large 
errors  of  our  analyses.  In  particular  the  height  analyses  did  not  fit  the  RAOBs 
closely  enough.  In  the  Northern  Hemisphere,  where  the  RAOBs  are  close 
together,  the  lack  of  fit  should  be  about  the  size  of  the  RAOB  observation  error 
standard  deviation  (OESD).  However  the  RAOB  OESD  at  500  hPa  is  only  about 
12  m,  while  our  typical  RAOB  analysis  difference  was  40  m. 

As  a  first  step  toward  addressing  this  problem,  we  removed  the  anchoring  of 
radiosonde  heights  to  the  first  guess  surface  pressure  and  conducted  a  test 
assimilation  for  24  hours.  Results  indicated  a  better  fit  to  the  500  hPa 
geopotential,  but  rapidly  growing,  large  errors  at  1000  hPa  as  a  result.  It 
became  clear,  then,  that  an  analysis  of  surface  pressure  was  required. 

A  new  surface  pressure  update  similar  to  that  used  by  the  ECMWF  was 
instituted.  The  new  formulation  uses  the  lowest  and  underground  layer 
height  increments  to  estimate  the  height  increment  at  o  =  1,  which  is  then 
converted  hydrostatically  to  a  pressure  increment.  With  these  corrections  the 
overall  impact  is  a  significant  improvement.  For  example.  Fig.  2  shows  the 
impact  on  the  500  hPa  height  analysis. 
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Figure  2:  Impact  of  improved  surface  pressure  analysis  on  500  hPa  rms 
analysis  error  for  the  Northern  Hemisphere  extratropics.  The  surface 
pressure  experiment  (dashed  line)  starts  on  Julian  day  325,  from  a  state 
identical  to  the  control  experiment  (solid).  The  RAOB  observation  error 
standard  deviation  (dotted)  is  shown  for  reference. 

3.1.13  Data  selection  for  bottom  layers 

Because  of  the  need  to  resolve  the  boundary  layer  structure  in  the  forecast 
model,  the  bottom  layers  of  the  model  (and  the  analysis)  are  very  thin.  This 
can  lead  to  problems  in  the  conversion  of  the  height  increments  to 
temperature  increments:  if  different  data  are  selected  for  adjacent  layers,  the 
height  increments  will  differ  by  an  amount  that  can  imply  an  unreasonably 
large  temperature  increment.  To  avoid  this  problem,  we  added  an  option  to 
the  OI  to  enforce  identical  data  selection  for  the  lowest  layers  (plus  the 
underground  layer).  This  was  done  by  temporarily  resetting  the  pressure  of 
the  sigma  layer  to  a  constant  value  for  all  these  layers  during  data  selection, 
and  subsequently  restoring  it  to  its  proper  value  when  the  actual  correlations 
for  the  analysis  weights  are  computed.  In  our  implementation,  this  data 
selection  was  enforced  for  the  lowest  4  layers. 

3.1.2  Radiance  retrieval 

A  major  disadvantage  of  the  operational  statistical  retrieval  is  that,  in  regions 
with  sufficient  conventional  data,  the  NWP  model  forecast  fields  themselves 
may  be  more  accurate  than  the  satellite  retrievals  for  initialization  purposes. 
Therefore,  assimilation  of  the  satellite  data  would  naturally  tend  to  degrade 
the  quality  of  forecasts  rather  than  to  improve  them.  This  property  of  the 
statistical  retrieval  is  due  to  its  dependence  on  precomputed  statistics  which 
imply  a  first  guess  which  is  essentially  climatology.  Isaacs'  unified  retrieval 
(Isaacs,  1987)  sought  to  circumvent  this  problem  by  using  the  statistical 
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retrieval  only  as  a  first  guess.  A  physical  retrieval  step  was  then  employed  to 
reconcile  the  retrieved  state  of  the  atmosphere  with  the  radiance  data  itself. 

Eyre  (1988)  points  out  several  other  problems  of  the  operational  retrieval 
schemes  in  addition  to  this  implied  first  guess  dependence.  These  include 
errors  from  the  cloud  clearing  algorithms  and  the  manner  in  which  the 
satellite  data  are  treated  in  the  assimilation  process.  For  CASCADE,  we  have 
used  Eyre's  optimal  retrieval  scheme,  which  operates  directly  on  the  radiance 
data,  thereby  bypassing  the  error  generating  preprocessing  and  cloud  clearing 
steps  in  the  conventional  approach.  His  method  may  be  considered  either  as  a 
retrieval  scheme  which  uses  the  forecast  model  to  provide  first  guesses  to  the 
desired  state  variables  and  then  provides  retrievals  which  are  consistent  with 
the  model,  or  alternatively,  and  equivalently,  as  a  module  within  the  analysis 
scheme  which  performs  the  analysis  of  the  satellite  data  in  the  vertical. 

The  background  or  first  guess  field  for  the  optimal  retrieval  of  Eyre  is  the 
forecast  model  itself.  In  our  implementation,  the  background  field  is  the 
analysis  produced  by  the  OI.  The  mathematics  of  the  retrieval  requires  the 
evaluation  of  the  forward  problem  and  its  derivatives  with  respect  to  the 
atmospheric  state  vector  elements.  The  optimal  retrieval  uses  information 
related  to  both  the  data  and  the  background  or  forecast  first  guess  field  to 
constrain  the  inversion. 


The  optimal  retrieval  is  characterized  mathematically  in  the  following  way: 
The  retrieval  seeks  to  find  the  most  probable  atmospheric  state  vector,  x, 
given  a  set  of  measurements,  i.e.  the  satel  ite  radiance  data,  y"».  By  Bayes’ 
theorem  this  conditional  probability,  is  proportional  to  the  product: 

P(y'"|x)-  P(x).  The  first  factor  is  the  probability  of  making  the  measurements 
given  the  state  of  the  atmosphere,  while  the  second  is  the  probability  of  the 
atmosphere  being  in  that  state.  The  former  can  be  related  to  the  covariance 
matrix  of  the  measurement  error  for  the  satellite  data,  while  the  latter  can  be 
related  to  the  forecast  model  error  covariance  assuming  that  the  NWP  model 
provides  the  first  guess.  It  is  assumed  that  both  error  sources  have  Gaussian 
error  statistics.  The  most  probable  atmospheric  state  given  the  measurements 
and  the  forecast  first  guess  then  involve  these  error  covariance  matrices  and, 
in  addition,  the  derivative  of  the  forward  problem  with  respect  to  the  state 
vectors. 


The  optimal  retrieval  scheme  treats  the  nonlinearities  of  the  retrieval 
problem  such  as  that  introduced  by  cloud  clearing  through  Newtonian 
iteration.  The  residuals,  i.e.  the  differences  between  the  actual  and  calculated 
radiances  for  a  given  guess  of  the  state  vector,  provide  information  on  the 
convergence  properties  of  the  retrieval.  Eyre  derives  a  retrieval  error 
covariance  matrix,  S,  which  can  be  used  to  address  the  problem  of  assigning 
weights  to  the  retrievals  within  the  horizontal  analysis.  The  problem  of 
correlation  between  errors  in  the  first  guess  (i.e.  the  forecast)  and  the 
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retrievals  is  also  addressed.  A  detailed  description  of  the  method  can  be  found 
in  Appendix  B. 

We  have  employed  the  optimal  retrieval  scheme  with  the  DMSP  microwave 
sensor  suite  of  instruments  (Falcone  and  Isaacs,  1987).  State  vector  elements 
which  can  be  treated  by  the  optimal  retrieval  include  temperature  and 
moisture  profiles,  surface  temperature,  and  cloud  parameters.  SSM/1  is 
particularly  useful  for  the  quality  control,  providing  a  precipitation  flag. 
Potentially,  rain  rate,  or  latent  heating  rate  information  from  SSM/I  could  be 
provided  to  the  initialization.  A  number  of  cloud  parameters  are  potentially 
retrievable  from  the  DMSP  sensors.  In  this  work  we  have  limited  the 
retrievals  to  temperature  and  moisture  profiles. 

In  the  configuration  provided  to  us  by  J.  Eyre  (TOVCFG),  the  scheme  was  set 
up  for  application  to  TIROS  Operational  Sunder  (TOYS)  data  and  included 
surface  microwave  emissivity,  cloud-top  pressure  and  effective  cloud  amount 
among  the  retrieved  parameters.  The  algorithm  used  a  forecast  profile  as  a 
background  constraint  and  as  a  first  guess  for  tempjerature  and  humidity 
parameters.  The  associated  error  covariance  matrix,  C,  were  constructed  from 
a  set  of  forecast  error  statistics  supplied  by  the  Forecasting  Research  Branch  of 
the  U.  K.  Meteorological  Office. 

The  TOVCFG  program  was  implemented  on  the  Alliant  computer  at  AER 
and  was  subsequently  modified  so  it  could  be  used  for  simulation  runs.  By 
being  able  to  directly  compare  the  retrieved  parameters  with  the  "true"  values 
used  to  generate  a  set  of  simulated  TOYS  radiances,  we  were  able  to  evaluate 
the  performance  of  the  retrieval  algorithm  and  at  the  same  time  check  for 
potential  coding  errors  which  might  have  occurred  in  the  process  of 
transferring  the  code  to  the  Alliant  system. 

As  part  of  our  simulation  toot,  we  developed  a  module  that  first  does  the 
vertical  interpolation  and  extrapolation  of  a  selected  forecast  profile  from  the 
nature  run  to  the  40  standard  National  Environmental,  Satellite,  and  Data 
Service  (NESDIS)  levels  used  by  the  radiative  transfer  algorithm.  A  "true" 
profile  X*  is  then  constructed  by  perturbing  the  resulting  background  profile 
as  follows: 

x' = X* + (1) 


where  A,  and  Vi  are  the  eigenvalues  and  associated  eigenvectors  of  the 
forecast  error  covariance  matrix  and  e,  are  samples  of  a  Gaussian  population 
with  zero  mean  and  unit  variance.  In  this  way,  the  perturbation  has  the  same 
covariance  as  the  forecast  errors.  Both  "true"  profile  and  background  profile 
are  checked  and  if  necessary,  corrected  for  supersaturation.  Finally,  a  vector  of 
simulated  TOYS  radiances  is  obtained  by  applying  the  forward  model  to  the 
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"true"  profile  and  by  adding  random  Gaussian  noise  of  covariance  E  to  the 
calculated  radiances. 

Figures  3  and  4  show  examples  of  simulated  retrievals  for  a  clear-sky  case  and 
a  cloudy  case,  respectively.  The  same  temperature  and  humidity  profiles  have 
been  used  in  both  experiments. 


Temperature  (K) 


Figure  3  Example  of  retrieval  for  a  clear  case.  The  solid  lines  are  the  first  guess 
errors  of  temperature  (left)  and  logarithm  of  specific  humidity  (right).  The 
dashed  lines  are  the  corresponding  retrieval  errors 
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Temperature  (K) 


Log  Mix.  Ratio  (g/kg) 


Figure  4  :  Same  as  Fig.  3,  but  for  a  cloudy  case. 


3.1.3  Horizontal  smoothing  algorithm 

A  horizontal  smoothing  algorithm  is  used  within  CASCADE  to  distribute 
retrievals  from  observing  locations  to  grid  points  in  the  horizontal.  We  have 
received  and  tested  the  recursive  filter  by  Hayden  and  Purser  (1988)  for  this 
purpose.  This  filter  is  a  fast  algorithm  that  is  essentially  equivalent  to  a 
Barnes  analysis.  We  generalized  the  software  to  allow  the  use  of  variable  grid 
spacing  (such  as  latitude-longitude  grids),  and  we  developed  guidelines  for 
choosing  the  adjustable  parameters  of  this  filter  based  on  the  work  of  Seaman 
(1989).  Seaman  derived  opHmal  values  for  the  Gaussian  length  scale  of  the 
Barnes  scheme  based  on  the  statistical  properties  of  the  background  field  and 
the  observations.  We  did  not,  however,  implement  this  part  of  CASCADE  in 
the  simulation  experiments  described  here. 

3.1.4  Qoud  data  processing 

The  purpose  of  this  module  is  to  use  cloud  data,  such  as  from  the  RTNEPH  or 
other  cloud  analysis  system,  to  refine  the  moisture  analysis.  It  is  also  used  to 
provide  latent  heating  rates  for  the  subsequent  initialization  step.  Both 
Molinari  (1982)  and  Fiorino  and  Warner  (1981)  obtained  improved  forecasts 
of  intensity  (relative  to  control  experiments)  when  specified  latent  heating 
was  used  during  a  12  h  assimilation  prior  to  the  forecast.  In  addition. 
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Molinari  obtained  improved  storm  track  forecasts  and  Fiorino  and  Warner 
obtained  improveo  precipitation  forecasts. 

Satellite-derived  cloud  parameters  have  great  potential  in  this  regard.  Several 
authors  (see  reviews  by  Martin  and  Scherer,  1973,  and  Atlas  and  Thiele,  1981) 
have  shown  that  precipitation  rates  can  be  specified  from  satellite  data,  and 
Julian  (1984)  has  proposed  that  the  vertical  profile  of  the  divergence  (or, 
equivalently,  vertical  velocity)  can  be  related  to  cloud  top  equivalent  black 
body  temperature  (EBBT)  in  the  tropics.  Therefore,  it  should  be  possible  to 
specify  latent  heating  rates  from  observed  cloud  characteristics.  In  our 
simulation  experiments  we  simply  used  the  precipitation  information 
existing  in  the  ECMWF  nature  run  as  a  proxy  for  cloud  data. 

A  possible  extension  of  this  work  would  be  to  use  cloudiness  data  in  another 
retrieval/analysis  step  in  the  framework  of  a  series  of  successive  optimal 
analyses/retrievals.  The  formalism  of  Eyre  (1988)  for  radiances  could  be 
applied  here  as  well,  with  the  forward  problem  consisting  of  the  computation 
of  cloud  value  from  the  analysis  variables. 

3.1.5  Initialization 

The  normal  mode  initialization  (NMI)  is  a  part  of  the  forecasting  system 
which  adjusts  the  initial  data  in  such  a  way  that  undesirable  gravity  waves  do 
not  grow  in  the  forecast.  The  ASAP  NMI  is  based  on  the  NMC  NMI,  which  is 
described  by  Ballish  (1980). 

If  diabatic  heating  effects  are  not  included,  NMI  tends  to  suppress  the  Hadley 
circulation.  A  number  of  devices  to  improve  the  initialized  tropical  fields 
within  the  context  of  adiabatic  NMI  have  been  suggested  (Puri,  1983;  Puri  et 
al.,  1982),  but  including  diabatic  processes  is  desirable,  because  the  latent  heat 
of  convection  is  so  clearly  important  in  the  balance  of  the  tropical  motions 
(Puri,  1985).  It  is  not  surprising,  therefore,  that  some  of  Rasch's  (1985b) 
findings  suggest  that  a  good  humidity  analysis  is  vital  to  diabatic  nonlinear 
NMI.  Unfortunately,  diabatic  nonlinear  NMI  has  been  found  to  have  poor 
convergence  characteristics  and  may  actually  increase  the  amplitudes  of  the 
high  frequency  noise  at  the  start  of  the  forecast  (Rasch,  1985a). 

An  additional  approach  to  improving  the  diabatic  nonlinear  NMI 
convergence  characteristics  is  to  hold  the  latent  heating  fixed.  Mohanty  et  al. 
(1986),  in  their  study  of  the  impact  of  diabatic  heating  on  nonlinear  NMI, 
obtained  their  best  results  using  "observed"  diabatic  heating.  These  were 
"observed"  by  calculating  the  residual  of  the  thermodynamic  equation 
averaged  over  24  h  from  ECMWF  analyses.  Similar  conclusions  were  drawn 
by  Puri  (1985). 

In  the  simulation  experiments  reported  here,  observed  precipitation  rates 
were  used  in  a  moist  initialization  procedure  which  consisted  of  an 
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initialization  step  and  a  forecast  modification  of  the  first  guess  forecast.  The 
forecast  modification  consisted  of  a  rescaling  of  the  predicted  convective 
heating  (and  precipitation)  by  the  ratio  of  observed  to  predicted  convective 
precipitation  for  all  grid  points  in  the  tropics  with  non  zero  predicted 
convective  rain.  The  factor  was  not  allowed  to  exceed  the  maximum  value  of 
10.  The  three  components  of  the  initialization  step  consisted  of  a  diabatic 
NMI,  a  separate  divergence  adjustment  step  in  which  vertical  profiles  of 
horizontal  divergence  were  adjusted  based  on  observed  precipitation  rates, 
and  a  temperature  and  moisture  adjustment  step  designed  to  bring  predicted 
convective  precipitation  rates  in  closer  agreements  with  observations.  The 
formulation  of  this  procedure  is  described  in  full  detail  in  Apjjendix  C. 

3.2  Experimental  set-up 

All  the  experiments  performed  in  this  project  were  done  in  simulation  mode. 
The  experimental  set-up  was  the  same  for  all  the  tests,  and  was  similar  to  the 
one  used  for  observing  system  simulation  experiments  (OSSEs)  in  our 
previous  work  (Grassotti  et  ai,  1989;  Hoffman  et  ai,  1990). 

The  data  is  extracted  from  a  nature  run,  which  is  an  extended  simulation 
performed  with  a  different  model.  From  an  initial  state  given  by  the  nature 
run  we  first  perform  what  we  call  a  spin-up  experiment.  It  consists  of  a  4  day 
forecast,  which  gives  a  simulated  atmospheric  state  with  errors  greater  than 
the  expected  analysis  errors,  followed  by  3  days  of  data  assimilation,  with  the 
baseline  system.  This  provides  the  initial  state  for  all  the  experiments.  Each 
experiment  consists  of  7  days  of  data  assimilation,  and  4  day  forecasts  starting 
from  days  3,  5  and  7  of  the  assimilation  period.  A  schematic  overview  of  this 
experimental  setup  is  shown  in  Fig.  5.  Since  the  "truth"  is  known,  as  given  by 
the  nature  run,  we  can  compute  analysis  and  forecast  errors,  by  subtracting 
the  nature  run  from  the  analyzed  and  forecast  fields. 

Forecasts 
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Initial  Spinup  Assimilation 

forecast  assimilation  experiment 


Fig.  5:  Schematic  showing  the  sequence  of  spinup  forecast,  spinup 
assimilation,  7-day  assimilation,  and  4-day  forecasts  from  day  3,  5, 
and  7. 
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3.2.1 


The  Nature  run 


The  nature  run  used  in  these  experiments  is  an  extended  simulation  created 
by  ECMWF.  Its  horizontal  resolution  corresponds  to  a  triangular  106 
truncation  in  sjjectral  space,  and  it  has  19  levels  in  the  vertical  (however,  the 
data  available  to  us  contained  data  on  14  mandatory  pressure  levels  between 
1000  hPa  and  10  hPa).  The  simulation  starts  from  OOGMT,  16  January  1979. 

The  tapes  containing  the  data  were  provided  to  us  by  NASA/GLA. 

In  addition  to  all  the  prognostic  variables  ^temperature,  winds,  specific 
humidity,  surface  pressure,  surface  temperature  and  moisture,  and  snow 
amount),  the  tapes  contain  a  number  of  diagnostic  fields  such  as  precipitation 
and  surface  fluxes.  The  data  are  stored  every  3  hours  for  31  days^ . 

3.2.2  The  simulated  observations 

The  simulated  observations  were  generated  by  adding  appropriate  errors  to 
"perfect  observations"  extracted  from  the  nature  run.  These  perfect 
observations  are  simply  the  value  of  the  nature  run  interpolated  to  the 
position  of  the  actual  observations  available  during  the  period  corresponding 
to  the  nature  run^.  We  did  not  interpolate  in  time.  Instead,  we  assigned  the 
observations  to  the  nearest  synoptic  time.  We  give  here  a  short  summary  of 
the  procedures  for  adding  errors,  which  are  described  in  detail  in  Appendices 
A  and  B. 

Errors  of  representativity  were  added  as  random,  uncorrelated  errors.  The 
error  magnitudes  were  determined  from  the  observed  spectra  reported  in 
Nastrom  and  Gage  (Nastrom  and  Gage,  1985;  Nastrom  et  al.,  1984).  No 
attempt  was  made  to  compensate  for  the  smoothing  of  the  highest  resolved 
wavenumbers  in  the  ECMWF  nature  run. 

We  added  random  errors  to  perfect  observations  of  types  1  (Raobs),  type  2 
(aircraft  obs),  and  type  6  (satellite  track  winds),  with  vertical  correlation  where 
appropriate.  For  type  4  (satellite  retrievals  from  TOYS),  we  performed  a  full 
simulation  rather  than  add  errors  to  the  perfect  observations.  The  error 
adding  software  additionally  sorted  the  observations  by  latitude  and,  in  the 
case  of  the  satellite  track  winds,  superobservations  were  formed. 

We  used  a  different  method  to  simulate  the  optimal  retrieval  of  SSM/T  data. 
The  SSM  retrievals  could  be  estimated  in  a  straightforward  fashion  by 
simulating  each  step  of  the  Eyre  procedure.  This  would  involve  computing 


^We  received  only  the  first  18  days  of  data. 

^The  perfect  observations  were  supplied  to  us  by  STC  . 
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radiances  from  the  nature  data  to  simulate  the  observations,  computing 
radiances  from  the  first  guess  (the  OI  analysis  in  our  case)  and  modifying  it 
iteratively  to  find  the  most  probable  profile  compatible  with  the  first  guess 
and  the  observations.  More  details  can  be  found  in  Appendix  B. 

An  accelerated  method  for  the  simulation  of  retrievals  of  profiles  and  other 
geophysical  parameters  from  microwave  sounders  has  been  used,  based  on 
Rodgers'  (1990)  linear  error  analysis.  The  approach  takes  advantage  of  the  fact 
that  the  target  profiles,  taken  from  the  nature  run  in  the  simulation  studies, 
are  known.  It  is  assumed  that  the  solution  of  the  non-linear  retrieval 
algorithm  converges  towards  a  point  close  to  the  target  profile  in  the 
"measurement  space"  and  that  the  magnitude  of  the  combined  radiometric 
and  forward  model  errors  remains  reasonably  small,  so  that  their  effect  can  be 
accounted  for  by  applying  a  linear  perturbation  to  the  target  profile.  Under 
these  assumptions,  the  solution  profile  can  be  expressed  as  a  linear 
combination  of  the  target  and  first  guess  profiles  and  the  radiometric  noise. 
The  weights  are  derived  from  the  error  covariance  matrices. 

The  amount  of  computation  required  to  simulate  a  retrieved  profile  with  this 
approach  is  roughly  equivalent  to  one  step  of  the  classical  Newton  iteration 
procedure  used  in  non-linear  retrieval  algorithms.  Computational  savings 
are  the  results  of  skipping  the  successive  iteration  steps  normally  needed  to 
bring  the  solution  to  convergence. 

An  assessment  of  the  validity  of  this  approach  was  made  based  on 
comparisons  between  the  solutions  obtained  with  the  accelerated  method  and 
the  ones  obtained  by  applying  Eyre’s  retrieval  algorithm.  The  Eyre  code  has 
been  modified  to  perform  retrievals  of  temperature  and  moisture  profiles 
from  simulated  SSM/1 ,  T  and  T-2  microwave  data.  Results  of  the 
comparison,  obtained  with  realistic  noise  data,  indicated  that  the  retrieval 
solution  can  be  estimated  using  the  fast  algorithm  with  an  accuracy  better 
than  1%. 

3.2.3  TTie  baseline  data  assimilation  system 

The  baseline  data  assimilation  system  consists  of  the  GSM,  ASAP  OI,  and  an 
adiabatic  NMI.  The  version  of  the  GSM  used  in  our  study  was  determined 
from  the  results  of  several  test  forecasts.  It  is  the  Phillips  Laboratory  GSM 
with  advanced  physics  (see  next  section).  We  used  it  in  a  configuration  with  a 
rhomboidal  40  truncation,  and  18  layers  (o  interface  values  are  1.0,  .990,  .973, 
.948,  .893,  .820,  .735,  .642,  .546,  .45,  .40,  .35,  .30,  .25,  .20,  .15,  .10,  .05,  and  0).  A 
time  step  of  15  minutes  was  used;  the  diffusion  coefficient  for  the  horizontal 
diffusion  was  set  at  2.55x10^-'’  m'l/s. 

The  ASAP  OI  system,  modified  as  described  above,  performed  the  analysis  on 
the  model  o  layers  on  a  81x82  Gau.ssian  grid. 


18 


The  adiabatic  NMI  was  identical  to  that  used  in  our  earlier  work.  It  is  a 
standard  2-iteration  Machenhauer  schen\e,  using  the  lowest  4  internal  modes 
with  a  frequency  cutoff  of  48  hours. 

3.2.4  The  spinup  experiment 

We  made  several  spin-up  forecasts,  with  different  versions  of  the  GSM.  We 
settled  on  a  version  that  uses  a  general  truncation  version  of  the 
hydrodynamics  code,  and  vectorized  versions  of  the  physics  packages 
(Nehrkorn  et  al.,  1991).  Its  physics  package  consists  of  the  Oregon  State 
University  planetary  boundary  layer  package  (Mahrt  et  al.,  1987),  the 
University  of  Illinois  Kuo  convection  (Soong  et  al.,  1985),  and  a  modification 
of  the  University  of  Utah  radiation  package  (Liou  et  al.,  1984),  developed  by 
the  PL  and  the  AER  climate  group.  It  includes  a  number  of  enhancements  to 
the  radiation  code  itself  and  uses  the  Slingo  (1987)  scheme  for  cloud  diagnosis, 
tuned  to  give  reasonable  cloud  cover  and  cloud  forcing  climatologies. 

Following  the  spinup  forecast,  a  three-day  assimilation  with  the  baseline  data 
assimilation  system  was  performed.  The  baseline  system  included  the  GSM, 
the  OI,  and  an  adiabatic  NMI,  as  described  above. 

3.3  The  CONTROL  OSSE 

This  OSSE  is  a  simple  7-day  continuation  of  the  spinup  assimilation 
experiment,  with  4-day  forecasts  generated  from  the  initialized  analyses  at 
days  3,  5,  and  7.  It  was  designed  to  document  the  performance  of  the  baseline 
data  assimilation  system,  and  to  provide  a  benchmark  for  comparison  with 
the  other  OSSEs. 

3.4  The  SSMTRAD  OSSE 

This  OSSE  differs  from  the  CONTROL  OSSE  only  in  that  the  radiance 
retrieval  step  of  SSM/T  radiances  was  included  in  the  data  assimilation 
scheme.  It  was  designed  to  test  the  radiaiice  retrieval  technique,  and  to 
determine  the  impact  of  the  SSM/T  data.  The  details  of  the  radiance  retrieval 
implementation  in  our  OSSE,  and  results  from  these  simulated  retrievals,  are 
discussed  in  detail  in  Appendix  B;  results  of  the  OSSE  itself  are  discussed  in 
section  4.2. 

3.5  The  MOIST  INIT  OSSE 

The  MOIST  INIT  OSSE  is  identical  to  CONTROL,  except  that  the  adiabatic 
NMI  is  replaced  by  a  moist  initialization  procedure.  The  moist  initialization 
was  designed  to  explore  and  test  methods  of  using  observed  precipitation 
rates  in  the  CASCADE  data  assimilation  scheme.  In  this  exploratory  study,  no 
errors  were  added  to  the  nature  run  precipitation.  The  impacts  thus  represent 
an  upper  limit  for  the  techniques  employed  here.  The  methods  were 
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implemented  as  an  initialization  step  and  a  modification  of  the  first  guess 
forecast. 

The  forecast  modification  was  motivated  by  the  finding  that  using  observed 
diabatic  heating  within  the  diabatic  NMI  lead  to  improved  results.  It  stands  to 
reason,  then,  that  additional  benefit  can  be  gained  by  forcing  the  forecast 
model  used  to  provide  the  first  guess  for  the  next  analysis  cycle  with  observed 
diabatic  heating  rates.  We  implemented  this  as  a  rescaling  of  the  predicted 
convective  heating  (and  precipitation),  at  each  time  step  and  at  all  grid  points 
in  the  tropics  with  non  zero  predicted  convective  rain,  by  the  ratio  of 
observed  to  predicted  convective  precipitation.  The  factor  was  not  allowed  to 
exceed  the  maximum  value  of  10.  The  observed  precipitation  was  the  nature 
run  total  precipitation  rate  for  the  6-hour  period  centered  at  the  inital  time 
(the  same  observed  precipitation  rate  was  used  for  the  forecast  modification 
and  the  moist  initialization).  In  a  real  data  experiment,  observed  precipitation 
rates  v/ould  have  to  be  inferred  from  mostly  remote  sensing  measurements. 

The  initialization  step  followed  all  other  analysis/retrieval  steps  of 
CASCADE.  The  three  components  of  the  initialization  step  consisted  of  a 
diabatic  NMI  (using  tendencies  computed  from  the  GSM  with  full  physics, 
with  diabatic  heating  in  the  tropics  rescaled  as  in  the  first  guess  forecast,  and  a 
2-iteration  Machenhauer  NMI),  a  separate  divergence  adjustment  step  in 
which  vertical  profiles  of  horizontal  divergence  were  adjusted  based  on 
observed  precipitation  rates,  and  a  temperature  and  moisture  adjustment  step 
designed  to  bring  predicted  convective  precipitation  rates  in  closer 
agreements  with  observations.  The  three  components  were  combined  into  a 
multistep  procedure,  as  follows: 

step  1:  diabatic  NMI 

step  2:  divergence  adjustment 

step  3:  temperature  and  moisture  adjustment 

step  4:  diabatic  NMI 

step  5:  temperature  and  moisture  adjustment 

The  development  and  testing  of  the  initialization  procedure  are  described  in 
detail  in  Appendix  C;  results  of  the  OSSE  itself  are  discussed  in  section  4.2. 

3.6  The  CASCADE  OSSE 

The  CASCADE  observing  system  simulation  experiment  combines  all  the 
parts  tested  separately  by  the  previous  two  OSSEs:  use  of  the  optimal  retrieval 
of  SSM/T  data  and  use  of  precipitation  data  both  for  modifying  the  moisture 
and  divergence  fields,  and  for  moist  initialization. 
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Results  of  Observing  System  Simularion  Experiments 
4.1  Performance  of  baseline  assimilation  system 

The  500  hPa  rms  errors  of  geopotential  height  (z)  displayed  in  Fig.  6  give  a 
good  illustration  of  the  spinup  phase  of  our  OSSEs.  During  the  initial  spinup 
forecast,  global  errors  grow  from  approximately  10  m  12  hours  into  the 
forecast  to  approximately  40  m  after  4  days.  The  spinup  assimilation,  which 
used  the  endpoint  of  the  4-day  forecast  as  its  starting  point,  reduces  the  error 
to  approximately  28  m  globally.  Analysis^  errors  have  stabilized  by  the  end  of 
the  3^ay  assimilation,  indicating  the  success  of  the  spinup  procedure.  In  the 
Northern  Hemisphere  extratropics,  where  data  density  and  quality  is  highest, 
the  error  reduction  is  quicker  and  final  rms  errors  are  on  the  order  of  18  m;  in 
contrast,  the  Southern  Hemisphere  shows  little  reduction  at  all  because  of  the 
lack  of  RAOB  data  (statistics  over  Australia,  not  shown  here,  show  a 
reduction  in  error  that  is  in  agreement  with  the  better  coverage  over  that 
region). 

The  analysis  errors  in  the  CONTROL  OSSE  (which  is  the  continuation  of  the 
spinup  assimilation  for  another  7  days)  show  little  variation  with  time,  aside 
from  a  slight  increase  in  Southern  Hemisphere  errors  in  the  second  half  of 
the  week. 

Time  series  of  700  hPa  relative  humidity  (rh)  and  the  700  hPa  and  300  hPa 
horizontal  wind  rms  errors  (Fig.  7-9)  show  qualitatively  similar  behavior  for 
CONTROL,  except  that  the  error  reduction  is  clearly  evident  only  in  the 
Northern  Hemisphere  extratropics.  In  the  tropics  and  the  Southern 
Hemisphere,  the  CONTROL  assimilation  merely  succeeds  in  halting  the 
growth  of  the  errors. 

Vertical  profiles  of  the  analysis  errors,  averaged  over  the  last  5  days  of  the 
assimilation,  provide  another  view  of  the  analysis  system  p>erformance.  The 
geopotential  errors  shown  in  Fig.  10  are  nearly  constant  with  height  below  the 
tropopause  with  values  between  20  m  (Northern  Hemisphere)  and  50  m 
(Southern  Hemisphere).  They  quickly  increase  with  height  in  the 
stratosphere;  this  increase  is  at  least  partially  due  to  errors  introduced  in  the 
interpolation  from  o  layers  to  mandatory  pressure  levels.  The  temperature 
errors  (Fig.  11)  show  a  similar  behavior  as  the  height  errors,  except  for 
comparatively  large  values  at  1000  hPa  and  850  hPa.  The  latter  are  most  likely 
due  to  different  treatment  of  underground  values  between  our  and  the 
ECMWF  postprocessor. 


^In  this  and  all  following  statistics,  analysis  refers  to  the  initialized  analysis. 
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6-  Root  mean  square  errors  of  500  hPa  geopotential  height  computed  for  the  globe  (a),  Northern  Hemisphere 
extratropics  (b),  Tropics  (c),  and  Southern  Hemisphere  extratropics  (d).  The  time  period  covers  the  sptnup 
forecast,  spinup  assimilation,  the  7-day  assimilation  and  the  4-day  forecasts  generated  from  it. 
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Fig.  7:  Root  mean  square  errors  of  700  hPa  relative  humidity  computed  for  the  globe  (a),  Northern  Hemisphere 
extratropics  (b),  Tropics  (c),  and  Southern  Hemisphere  extratropics  (d).  The  time  period  covers  the  spinup 
forecast,  spinup  assimilation,  the  7-day  assimilation  and  the  4-day  forecasts  generated  from  it. 
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Fig.  9:  Root  mean  square  errors  of  300  hPa  horizontal  wind  computed  for  the  globe  (a),  Northern  Hemisphere 
extratropics  (b),  Tropics  (c),  and  Southern  Hemisphere  extratropics  (d).  The  time  period  covers  the  spinup 
forecast,  spinup  assimilation,  the  7-day  assimilation  and  the  4-day  forecasts  generated  from  it. 
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Fig.  10:  Vertical  profiles  of  the  root  mean  square  errors  of  geopotential  height  averaged  over  the  last  5  days  of  the 
assimilation,  for  the  globe  (a),  Northern  Hemisphere  extratropics  (b),  Tropics  (c),  and  Southern  Hemisphere 
extratropics  (d). 
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Fig.  11:  Vertical  profiles  of  the  root  mean  square  errors  of  temperature  averaged  over  the  last  5  days  of  the  assimilation, 
for  the  globe  (a),  Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and  Southern  Hemisphere  extratropics  (d). 


Relative  humidity  errors  (Fig.  12)  are  smallest  near  the  surface,  and  show 
little  variation  between  700  hPa  and  300  hPa.  Rh  values  above  300  hPa  were 
available  both  from  the  nature  run  and  our  analyses,  but  since  there  are 
virtually  no  data  at  those  levels  those  values  are  not  displayed  here. 

The  error  profile  for  the  horizontal  wind  (Fig.  13)  has  a  distinct  maximum  at 
the  jet  level,  which  is  most  pronounced  in  the  extratropics.  Error  magnitudes 
at  that  level  are  approximately  7  m/s  in  the  Northern  Hemisphere,  and  11 
m/s  in  the  Southern  Hemisphere. 

The  performance  of  the  CONTROL  data  assimilation  system  can  be  compared 
with  that  of  previous  OSSEs  (Hoffman  et  ai,  1989;  Hoffman  et  al.,  1990, 
Grassotti  et  ai,  1989;  and  Grassotti  et  ai,  1991,  denoted  H89,  H90,  G89,  and  G91 
in  the  following)  using  the  Phillips  Laboratory  data  assimilation  system.  The 
first  set  of  these  OSSEs,  denoted  STATSAT,  SSM,  SSM+TOVS,  WINDSAT  in 
H89,  G89,  and  H90,  was  designed  to  assess  the  impact  of  statistical  retrievals  of 
temperature  and  moisture  from  SSM/T  measurements,  both  in  combination 
with  and  as  a  replacement  for  TOYS  statistical  retrievals,  and  of  an  active 
Doppler  Lidar  wind  sounder.  The  second  set,  denoted  CONTROL  and  LIDAR 
in  G91,  used  the  same  Ol,  but  with  the  same  improvements  in  the  use  of 
RAOB  data  and  the  surface  pressure  analysis  that  are  being  utilized  in  our 
present  OSSEs  (the  LIDAR  OSSE  used  simulated  data  from  a  scaled  down 
version  of  the  instrument  simulated  in  WINDSAT).  Neither  of  the  two 
previous  sets  used  the  new  physics  packages  in  the  GSM,  and  they  both  were 
run  at  an  R30  and  12-layer  resolution. 

The  present  CONTROL  OSSE  corresponds  most  closely  to  the  STATSAT 
experiment  of  G89  and  the  CONTROL  experiment  of  G91 .  The  results  of  the 
two  CONTROL  experiments  are  quite  similar:  500  hPa  height  and  850  hPa  and 
200  hPa  zonal  wind  analysis  errors  show  identical  behavior  and  have  roughly 
the  same  magnitudes.  Rh  analysis  errors  are  slightly  smaller  in  the  present 
CONTROL,  particularly  in  the  tropics,  because  of  the  improved  forecast 
model  (since  only  RAOB  data  is  used  in  the  moisture  analysis,  the  first  guess 
field  has  a  large  effect  on  the  analysis).  By  contrast,  the  STATSAT  experiment 
described  in  H89  suffered  from  a  suboptimal  use  of  RAOB  data  and  a  poor 
surface  pressure  analysis,  resulting  in  large  analysis  errors:  height  errors  of  35 
m  at  500  hPa,  compared  to  28  m  in  the  present  CONTROL;  zonal  wind  errors 
approximately  Im/s  larger  than  in  CONTROL;  and  rh  errors  approximately 
2%  larger  than  in  CONTROL. 
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Fig.  12:  Vertical  profiles  of  the  root  mean  square  errors  of  relative  humidity  averaged  over  the  last  5  days  of  the 
assimilation,  for  the  globe  (a),  Northern  Hemisphere  ext  rat  ropics  (b),  Tropics  (c),  and  Southern  Hemisphere 
extratropics  (d). 
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Fig.  13:  Vertical  profiles  of  the  root  mean  square  errors  of  the  horizontal  wind  (vector  difference)  averaged  over  the  last 
5  days  of  the  assimilation,  for  the  globe  (a),  Northern  Hemisphere  extratropics  (b),  Tropics  (c),  and  Southern 
Hemisphere  extratropics  (d). 


4.2  OSSE  impacts 


In  the  following  discussion  of  analysis  errors,  one  must  bear  in  mind  that 
these  analyses  have  been  initialized,  either  by  the  adiabatic  NMI  (for  the 
CONTROL  and  SSMTRAD  OSSEs),  or  by  the  moist  initialization  procedure 
(for  the  MOIST  INIT  and  CASCADE  OSSEs).  The  impacts  seen  in  the  analysis 
errors  are  thus  due  to  the  combined  effects  of  the  improved  first  guess  and 
additional  data  available  to  the  analysis,  and  the  initialization  procedure.  The 
regions  for  which  statistics  are  presented  in  the  following  are  defined  such 
that  "Tropics"  only  includes  points  between  20°S  and  20°N;  since  the  moist 
initializ.ation  procedure  Wcis  applied  to  points  equatorward  of  30°,  this  also 
affected  statistics  for  the  regions  denoted  "extratropics”,  which  include  points 
between  20°  and  80°. 

The  addition  of  SSM/T  data  (in  the  SSMTRAD  OSSE)  and/or  the  moist 
initialization  procedure  (in  the  MOIST  INIT  and  CASCADE  OSSEs)  have  no 
effect  on  the  500  hPa  height  analysis  errors  (Fig.  6).  From  Fig.  10  it  can  be  seen 
that  this  is  true  for  all  levels  at  and  below  500  hPa;  outside  the  Northern 
Hemisphere,  there  is  a  slight  positive  impact  of  the  SSM/T  data  below  100 
hPa.  However,  a  substantial  negative  impact  is  evident  above  100  hPa.  The 
corresponding  profiles  of  temperature  (Fig.  11)  show  zero  global  impact  of  the 
moist  initialization  procedure,  with  a  slight  improvement  in  the  tropics 
below  300  hPa,  and  a  slight  degradation  in  the  Southern  Hemisphere.  The 
SSM/T  data  has  little  or  no  effect  below  300  hPa,  a  positive  impact  at  100  and 
30  hPa,  but  a  small  negative  impact  at  250  hPa,  and  a  substantial  negative 
impact  at  70  and  50  hPa.  These  SSM/T  impacts  are  consistent  with  the 
retrieval  error  statistics  discussed  in  Appendix  B,  except  that  the  degradation 
in  the  stratosphere  introduced  at  each  retrieval  time  period  is  smaller  than 
the  accumulated  effect  evident  in  the  SSMTRAD  analysis  errors. 

The  relative  humidity  errors  show  a  clear  and  positive  effect  of  the  additional 
SSM/T  data,  both  in  the  700  hPa  time  series  (Fig.  7)  and  in  the  vertical  profiles 
of  rh  error  (Fig.  12).  Rh  errors  are  reduced  by  up  to  10%  (from  30%  to  20%)  in 
the  Southern  Hemisphere,  whereas  impacts  are  smaller  (around  5%)  in  the 
Northern  Hemisphere.  Interestingly,  the  rh  error  reduction  computed  from 
the  retrieval  statistics  (see  Appendix  B)  shows  reductions  that  increase  with 
height,  whereas  the  cumulative  effects  shown  in  Fig.  12  are  largest  at  700  and 
500  hPa.4 


^The  rh  analyses  in  the  stratosphere  (not  shown  here)  arc  generally  poorer  in  the  SSMTRAD 
OSSE.  Since  the  SSM/T  retrievals  of  moisture  arc  only  available  at  and  below  300  hPa,  this 
degradation  is  caused  by  changes  to  the  temperature  and,  indirectly,  the  wind  fields. 
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The  moist  initialization  has  a  very  small  impact  on  rh  errors  in  the  tropics  in 
the  MOIST  INIT  OSSE  (less  than  a  1%  error  reduction).  From  Fig.  7  it  can  be 
seen  that  a  larger  impact  exists  during  the  first  3  days  of  MOIST  INIT.  In  the 
CASCADE  OSSE,  the  positive  impact  persists  throughout  the  seven  day 
period,  and  the  CASCADE  OSSE  consequently  has  the  lowest  rh  errors  in  the 
tropics.  Outside  the  tropics,  there  is  little  difference  between  the  SSMTRAD 
and  CASCADE  rh  errors. 

The  largest  effect  of  the  moist  initialization  is  evident  in  the  tropical  wind 
fields  (Fig.  13),  where  errors  are  reduced  by  as  much  as  2  m/s  (MOIST  INIT, 
CASCADE  vs.  CONTROL,  SSMTRAD).  Improvements  are  smaller,  but  still 
noticeable,  in  the  extratropics.  This  is  consistent  with  the  off-line  tests  of  the 
moist  initialization  procedure  discussed  in  Appendix  C.  In  the  Southern 
Hemisphere,  the  SSM/T  data  also  improves  the  wind  field  below  150  hPa 
(SSMTRAD  vs.  CONTROL).  A  degradation  of  the  horizontal  winds  is  evident 
in  SSMTRAD  at  and  above  the  tropopause  in  all  regions;  this  is  related  to  the 
degradation  of  the  temperature  and  height  fields  at  those  levels  discussed 
above.  The  wind  errors  of  the  CASCADE  OSSE  are  in  general  equal  to  the 
smaller  of  SSMTRAD  or  MOIST  INIT,  except  in  the  Northern  Hemisphere 
where  CASCADE  errors  are  slightly  smaller  than  any  of  the  other  OSSEs. 

Time  series  of  the  700  hPa  and  300  hPa  horizontal  wind  errors  (Fig.  8  and  9) 
show  that  effects  evident  in  the  5-day  averages  are  typical  of  the  impacts 
throughout  the  7-day  assimilation. 

In  summary,  the  analysis  errors  indicate  a  small  positive  impact  of  SSM/T 
data  on  tropospheric  temperatures  in  the  tropics  and  Southern  Hemisphere,  a 
positive  impact  on  trop>ospheric  winds  in  the  Southern  Hemisphere,  and  a 
large  positive  impact  on  rh  in  the  tropics  and  the  Southern  Hemisphere.  A 
negative  impact  of  the  SSM/T  data  is  evident  in  the  temperature,  height,  and 
wind  fields  in  the  lower  stratosphere.  The  moist  initialization  lowered  errors 
in  the  wind  fields,  particularly  in  the  tropics,  and  rh  errors  in  the  tropics. 

With  the  exception  of  the  lower  stratosphere,  the  CASCADE  OSSE  resulted  in 
the  lowest  analysis  errors  overall,  both  because  positive  impacts  from  the 
.MOIST  INIT  and  SSMTRAD  complement  each  other,  and  because  in  some 
cases  (e.g.,  tropical  rh)  positive  impacts  are  partially  additive. 

The  analysis  errors  discussed  up  to  now,  while  of  interest  by  themselves,  also 
affect  the  quality  of  the  ensuing  forecasts.  The  forecast  errors  for  the 
individual  forecasts  are  shown  in  Fig.  6-9,  and  average  errors  from  the  day  3, 
5,  and  7  forecasts  are  shown  in  Fig.  14-  17.  As  might  be  expected  from  the 
analysis  error  impacts,  there  is  no  consistent  impact  on  the  500  hPa  height 
errors:  while  the  averaged  forecast  errors  show  a  positive  impact  for 
SSMTRAD  and  MOIST  INIT,  examination  of  the  individual  forecasts 
indicates  near  zero  or  negative  impacts  for  one  or  two  of  the  three  forecasts. 

Rh  errors  at  850  hPa  show  positive  impacts  of  SSM/T  data  in  the  Southern 
Hemisphere,  where  SSMTRAD  and  CASCADE  are  smaller  by  5%  (at  the 
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initial  time)  to  1%  (at  3  days).  By  the  fourth  day  of  the  forecast  there  is  no 
impact.  In  the  tropics,  CASCADE  forecast  errors  are  smaller  than  either 
SSMTRAD  or  MOIST  INIT  during  the  first  1  to  2  days  of  the  forecast.  The  rh 
forecast  improvements  are  evident  in  all  three  forecasts. 

Finally,  forecast  errors  of  zonal  wind  at  850  hPa  and  200  hPa  reflect  the 
improvements  to  the  Southern  Hemisphere  wind  analyses  in  SSMTRAD, 
and  the  tropical  wind  analyses  in  MOIST  INIT.  The  latter  does  not  persist  past 
the  first  day  of  the  forecast,  whereas  the  Southern  Hemisphere  850  hPa  wind 
field  errors  are  smaller  throughout  the  forecast  period  in  SSMTRAD.  In  the 
Northern  Hemisphere  the  averaged  forecast  errors  seem  to  indicate  a  positive 
impact  beyond  2  days,  but,  as  in  the  case  of  the  500  hPa  heights,  this  impact  is 
not  evident  in  all  3  forecasts  and  is  most  likely  due  to  sampling  limitations. 

Compared  to  our  previous  OSSEs  using  SSM/T  data,  impacts  seeii  here  are 
generally  smaller  due  to  the  smaller  errors  of  the  CONTROL  experiment.  In 
particular,  the  small  positive  impact  on  Northern  Hemisphere  heights  (2-3 
m)  seen  in  G89  is  absent  here,  and  rh  analysis  error  improvements  of  5%  are 
only  3-4%  here;  forecast  rh  error  improvements  in  G89  decreased  to  2%  at  day 
4,  compared  to  0%  at  day  2  here;  however,  the  absolute  error  magnitudes  are 
smaller  because  of  the  improved  analysis  and  forecast  model.  Rh  forecast 
errors  in  particular  were  saturating  near  30%  in  G89,  whereas  the  present 
error  curves  have  not  flattened  out  by  day  4,  when  they  reach  28%. 
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Fig.  14:  Root  mean  square  errors  of  500  hPa  geopotential  height  as  a  function  of  forecast  lead  time,  averaged  over  the 
day  3, 5,  and  7  forecasts,  for  the  globe  (a).  Northern  Hemisphere  extratropics  (b).  Tropics  (c),  and  Southern 
Hemisphere  extratropics  (d). 
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Fig.  15:  Root  mean  square  errors  of  850  hPa  relative  humidity  as  a  function  of  forecast  lead  time,  averaged  over  the  day 
3, 5,  and  7  forecasts,  for  the  globe  (a).  Northern  Hemisphere  extratropics  (b),  Tropics  (c),  and  Southern 
Hemisphere  extratropics  (d). 
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Fig.  16:  Root  mean  square  errors  of  850  hPa  zonal  wind  as  a  function  of  forecast  lead  time,  averaged  over  the  day  3, 5, 
and  7  forecasts,  for  the  globe  (a),  Northern  Hemisphere  extratropics  (b),  Tropics  (c),  and  Southern  Hemisphere 
extratropics  id). 


Global  rms  200  mb  error  (averaged  lest  error) 


5  Summary  and  Conclusions 

A  set  of  OSSEs  was  performed  to  assess  the  impact  of  SSM/T  temperature  and 
moisture  data,  and  of  a  moist  initialization  procedure  making  use  of  observed 
precipitation  rates.  The  SSM/T  data,  as  well  as  all  the  conventional  data,  were 
simulated  with  as  realistic  errors  as  possible.  The  observed  precipitation  rates, 
on  the  other  hand,  did  not  include  any  observational  errors,  and  the  results 
thus  represent  an  upper  limit  of  the  positive  impact  of  our  particular 
initialization  procedure. 

The  analysis  system  used  in  the  CONTROL  experiment,  which  used  the  full 
complement  of  conventional  data  in  an  OI  assimilation  scheme,  performed 
satisfactorily.  Comparisons  with  similar  assimilation  experiments  using 
earlier  versions  of  the  data  assimilation  system  showed  the  beneficial  effects 
of  improvements  to  the  Ol  (better  use  of  RAOB  data  and  improved  surface 
analysis)  and  the  GSM  (higher  resolution  and  more  sophisticated  physics). 

The  analysis  errors  indicate  a  small  positive  impact  of  SSM/T  data  on 
tropospheric  temperatures  in  the  tropics  and  Southern  Hemisphere,  a 
positive  impact  on  tropospheric  winds  in  the  Southern  Hemisphere,  and  a 
large  positive  impact  on  rh  in  the  tropics  and  the  Southern  Hemisphere.  A 
negative  impact  of  the  SSM/T  data  is  evident  in  the  temperature,  height,  and 
wind  fields  in  the  lower  stratosphere. 

The  moist  initialization  lowered  analysis  errors  in  the  wind  fields, 
particularly  in  the  tropics,  and  rh  errors  in  the  tropics.  With  the  exception  of 
the  lower  stratosphere,  the  CASCADE  OSSE,  which  combined  the  SSM/T  data 
and  the  moist  initialization,  resulted  in  the  lowest  analysis  errors  overall, 
both  because  positive  impacts  from  the  MOIST  INIT  and  SSMTRAD 
complement  each  other,  and  because  in  some  cases  (e.g.,  tropical  rh)  positive 
impacts  are  partially  additive. 

Forecast  errors  showed  positive  impacts  in  only  selected  variables  and 
geographical  regions.  Rh  errors  at  850  hPa  show  positive  impacts  of  SSM/T 
data  in  the  Southern  Hemisphere,  where  SSMTRAD  and  CASCADE  are 
smaller  by  5%  (at  the  initial  time)  to  1%  (at  3  days).  By  the  4‘h  day  of  the 
forecast  there  is  no  impact.  In  the  tropics,  CASCADE  forecast  errors  are 
smaller  than  either  SSMTRAD  or  MOIST  INIT  during  the  first  1  to  2  days  of 
the  forecast.  Finally,  forecast  errors  of  zonal  wind  at  850  hPa  and  200  hPa 
reflect  the  improvements  to  the  Southern  Hemisphere  wind  analyses  in 
SSMTRAD,  and  the  tropical  wind  analyses  in  MOIST  INIT.  The  latter  does 
not  persist  past  the  first  day  of  the  forecast,  whereas  the  Southern  Hemisphere 
850  hPa  wind  field  errors  are  smaller  throughout  the  forecast  period  in 
SSMTRAD. 
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Compared  to  our  previous  OSSEs  using  SSM/T  data,  impacts  seen  here  are 
generally  smaller  due  to  the  smaller  errors  of  the  CONTROL  experiment.  In 
particular,  the  small  positive  impact  on  Northern  Hemisphere  heights  (2-3 
m)  seen  in  G89  is  absent  here,  and  rh  analysis  and  forecast  error 
improvements  are  smaller  here;  however,  the  absolute  error  magnitudes  are 
smaller  because  of  the  improved  analysis  and  forecast  model.  The  negative 
impact  of  SSM/T  data  in  the  stratosphere  was  absent  in  our  previous  OSSE. 

It  is  remarkable  that  even  the  relatively  crude  moist  initialization  procedure 
implemented  in  the  data  assimilation  scheme  in  MOIST  EMIT  and  CASCADE 
results  in  the  analysis  and  initial  forecast  error  improvements  seen  here, 
particularly  since  the  spinup  behavior  of  precipitation  shown  in  Appendix  C 
is  far  from  perfect.  This  suggests  that  additional  improvements  might  be 
possible  if  more  sophisticated  procedures  are  employed,  and  similar 
approaches  applied  to  large-scale  precipitation  outside  the  tropics.  It  remains 
to  be  determined,  however,  how  much  the  impacts  will  be  modified  by  the 
inclusion  of  realistic  measurement  errors. 
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Appendix  A:  Data  Simularion  Methodology  for  the  Standard 
Observations 

A.l  Representativeness  of  the  nature  data. 

Scales  on  the  order  of  15  minutes  to  an  hour  or  two  are  not  in  the  nature  data 
because  of  model  resolution;  similarly  for  scales  from  25  km  to  100  km.  These 
scales  do  affect  the  observing  systems  (RAOBs  and  SATEMs  for  example). 
These  scales  should  be  "added"  to  the  nature  data  before  sampling  it  or 
otherwise  simulated.  Simulating  this  effect  essentially  adds  errors  to  the 
observations.  We  have  added  a  random  component  of  the  proper  size  to  each 
perfect  observation.  These  need  not  be  correlated  horizontally  since  the 
correlation  length  scale  is  so  small  for  these  scales.  It  is  also  unclear  whether 
these  errors  should  have  any  vertical  correlations,  so  we  have  not  introduced 
any  vertical  correlations.  The  proper  size  of  these  errors  was  determined 
using  an  energy  spectrum  of  E=de/dk=c  klf,  with  y=-5/3,  deduced  from 
observed  atmospheric  spectra  (Nastrom  and  Gage,  1985;  Nastrom  et  al.,  1984). 
Since  different  instruments  average  over  different  space-time  volumes, 
differently  sized  errors  were  used  for  each  data  type. 

The  proportionality  constant  c  was  obtained  by  fitting  this  equation  to  values 
in  Nastrom's  table,  using  entries  for  all  flights.  The  variance  under  this  curve 
was  then  determined  between  the  wavenumbers  2n/L\  and  In/Li,  where  Li  is 
the  resolution  of  the  nature  run  (T106  =  310  km),  and  L2  the  averaging  length 
scale  of  the  instrument.  There  is  less  than  a  10%  change  in  the  standard 
deviation  for  L2  going  from  20  km  to  0  km.  Resulting  values  are: 

winds  =  1.9  m/s  for  troposphere  (below  100  hPa),  2.0  m/s  (100  hPa  and 
above) 

temperature  =  .8  K  for  troposphere  ,  1.8  K  for  stratosphere  (this  value  was 
not  used,  because  it  might  have  been  caused  by  aircraft  bobbing) 

No  attempt  was  made  to  compensate  for  the  smoothing  of  the  highest 
resolved  wavenumbers  in  the  ECMWF  nature  run.  It  was  decided  that  the 
uncertainty  in  the  observed  spectra  did  not  warrant  this.  The  temperature 
values  are  not  used  explicitly  in  the  simulation  of  the  RAOBs,  although  the 
temperature  errors  resulting  from  the  specified  height  errors  seem  to  be  of  the 
right  magnitude. 

A.2  Time  interpolation. 

The  perfect  data  that  are  being  sampled  are  at  the  synoptic  times  (0, 6, 12  and 
18  GMT)  for  the  data  within  plus  or  minus  3  hours  of  these  times.  To 
interpolate  properly  in  time  one  would  need  to  do  twice  as  much  work  as 
without  time  interpolation.  For  example  for  AIREPs  at  4  GMT  one  would 
need  to  create  perfect  data  at  0  GMT  for  the  6  GMT  data  locations.  Similarly 
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for  reports  at  8  GMT  we  would  also  need  to  process  the  6  GMT  data  against 
the  12  GMT  nature  file.  It  might  be  sufficient  to  interpolate  the  12  -  00  GMT 
nature  run  difference  to  the  6  GMT  observations. 

It  might  be  best  to  perform  such  a  time  interpolation.  It  would  also  be  better  to 
improve  ASAP  to  calculate  observation  increments  which  interpolate  the 
first  guess  to  each  observing  time.  This  is  the  nominal  but  unaffordable  case, 
and  there  are  several  alternative  approaches.  The  first  alternative  is  to  do  just 
the  time  interpolation  and  not  modify  ASAP.  A  second  alternative  is  to 
simulate  the  effect  of  the  time  interpolation  by  adding  errors  to  the 
observations  of  the  same  size  as  those  introduced  by  ASAP  in  not  performing 
its  own  time  interpolation.  ASAP  would  also  be  left  unchanged  in  this 
scenario.  The  third  alternative  is  to  simulate  improving  ASAP  by  not  adding 
any  time  interpolation  errors. 

If  we  followed  the  second  alternative,  we  could,  as  an  approximation,  increase 
the  magnitude  of  the  errors  associated  with  for  the  off-time  observations.  Of 
course  the  errors  due  to  timeliness  are  correlated.  The  size  of  these  errors 
could  be  estimated  by  processing  one  time  period  including  time 
interpolation.  Alternatively,  we  could  examine  the  size  of  persistence  forecast 
errors  for  1,  2  and  3  hours.  This  would  be  simpler  and  give  us  stable  statistics 
at  all  latitudes  and  over  the  ocean.  The  original  simulated  errors  could  then 
be  amplified  by  different  factors  for  different  time  lags  and  latitude  (differently 
for  ocean  and  land  perhaps)  so  that  the  standard  deviations  are  correct. 

After  discussion  with  the  contract  monitor,  we  decided  to  adopt  the  third 
alternative,  which  is  by  far  the  simplest  to  implement. 

A.3  Measurement  error 

The  measurement  error  varies  with  data  type,  and  they  are  simulated  as 
follows: 

A.3.1  Type  1:  Raobs,  Pibal,  etc. 

Errors  are  expected  to  be  uncorrelated  horizontally  between  pairs  but 
correlations  in  the  vertical  are  expected.  We  obtained  the  correlations  from 
the  statistics  used  in  the  ECMWF  OI  scheme.  We  performed  an  eigenvector 
(or  EOF)  analysis  of  the  covariance  matrix  to  simulate  the  height  errors.  From 
those,  layer  temp)erature  errors  were  derived  hydrostatic^  and  mandatory 
and  significant  level  temperature  error  were  derived  that  were  consistent 
with  the  layer  temperature  errors.  Wind  errors  were  specified  using  the  EOFs 
of  its  vertical  covariance  matrix,  and  interpolated  to  significant  levels  as 
needed.  Rh  errors  were  specified  as  uncorrelated  random  numbers,  using 
instrument  errors  from  Elliott  and  Gaffen  (1991)  plus  an  assumed  10%  error 
of  representativity.  Note  that  these  instruments  do  not  sample  nature 
instantaneously.  A  radiosonde  ascent  to  100  hPa  takes  on  the  order  of  one 
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hour  or  more.  However,  as  discussed  above,  we  did  not  explicitly  treat  the 
time  interpolation  errors  of  any  of  the  observations. 

A.3.2  Type  2:  Aireps. 

Potential  height  assignment  errors  and  instrument  biases  can  introduce 
horizontally  correlated  errors  for  individual  aircraft.  However,  aircraft  ID  is 
not  saved  in  our  data  sets  so  it  would  be  difficult  to  simulate  this.  Since  the 
numbers  of  observations  is  so  small,  and  there  is  also  little  guidance  in  the 
choice  of  proper  correlation  lengths,  we  just  used  random  Gaussian  errors  for 
u  and  V  and  T.  Instrument  errors  were  taken  from  Hoffman  et  al.  (1990) 

A.3.3  Type  4:  Satems. 

Statistical  satellite  retrievals  for  TOVS  were  calculated  in  simulation  by 
performing  a  forward  problem  computation  and  a  statistical  retrieval. 
Retrievals  were  only  computed  over  water,  and  equatorward  of  65  degrees  (to 
avoid  complications  of  surface  emissivity  over  land  and  ice). 

The  forward  model  code  was  obtained  from  jon  Eyre.  The  input  profile  to  the 
forward  problem  consists  of  nature  run  temperature  and  rh  profiles  (at  the  14 
mandatory  levels),  bilinearly  interpolated  to  the  sounding  location  from  a  1 
by  1  degree  nature  run  data  set  and  surface  values  of  pressure,  temperature, 
rh,  and  u  and  v  (the  latter  obtained  using  routine  GETSFC  from  the  GSM 
preprocessor),  along  with  an  effective  cloud  cover.  The  cloud  cover  is  taken  as 
the  maximum  of  the  parameters  PCCl  and  PCC2  from  the  original  (real) 
retrievals.  Errors  of  represen tativity  appropriate  for  the  125  km  averaging 
length  scale  of  the  retrievals  are  added  to  the  nature  run  data  before  input  to 
the  forward  problem.  The  mandatory  level  data  is  then  interpolated  and 
extrapolated  to  the  40  (15  for  moisture)  forward  problem  levels.  The 
interpolation  matrices  have  been  derived  from  a  data  base  of  retrieval  profiles 
provided  by  McMillan  (NESDIS).  Instrument  errors  are  added  to  the 
computed  radiances. 

The  statistical  retrieval  used  a  D-matrix  that  has  been  computed  from  an 
independent  sample  of  nature  run  data  (the  4  days  from  Jan.  16  to  Jan.  20, 

1979)  Separate  D-matrices  were  derived  for  "clear"  and  "cloudy"  retrievals, 
but  no  further  stratification  was  done.  The  geometry  is  not  explicitly 
simulated:  the  forward  problem  was  done  for  the  HIRS  and  MSU  channels, 
using  the  exact  same  profile,  and  nadir  conditions  were  assumed  for  all 
retrievals.  Separation  between  clear  and  cloudy  retrievals  were  based  on  a 
single  cutoff  value  for  cloud  cover,  and  not  the  more  complicated  criterion 
used  in  the  actual  retrievals  (based  on  the  number  of  cloud-free  HIRS  IFOVs); 
similarly,  the  "mixed"  retrieval  were  not  simulated.  The  threshold  was 
selected  based  on  the  reported  distribution  of  cloudy  and  clear  retrievals,  and 
the  distribution  of  max(PCCl,  PCC2). 
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The  retrieved  profiles  of  temp)erature  were  then  converted  to  layer  mean 
temperatures  and/or  layer  thicknesses,  depending  on  the  format  of  the  perfect 
observation,  and  layer  precipitable  water.  Note  that  the  moisture  retrievals 
are  not  used  by  our  OI. 

A.3.4  Type  6:  CDWs. 

For  satellite  winds,  horizontally  filtered  random  instrument  errors  were 
added,  as  detailed  below.  Error  magnitudes  were  taken  from  Hoffman  et  al. 
(1990).  Superobservations  were  formed. 

The  CDWs  simulated  by  Dey  et  al.  (1985)  have  random  error  only,  whereas  we 
know  a  sizable  fraction  of  the  error  is  due  to  height  assignment  errors 
(McPherson,  1984).  Further  sources  of  error  are  the  fact  that  cloud  base  not  the 
cloud  top  may  be  the  best  level  for  assigning  the  winds  and  the  deviations 
from  unit  emissivity  for  the  cloud.  These  errors,  especially  the  height  errors 
tend  to  be  very  well  correlated  (at  least  at  a  particular  synoptic  time)  for  a 
particular  data  producer.  The  cause  of  the  height  assignment  errors  is 
fundamental.  Generally  two  approaches  to  height  assignment  have  been 
used: 

1  The  climatological  approach.  For  example,  assigning  all  low  level  wind 
to  900  hPa  or  assigning  all  high  level  wind  to  the  climatological 
tropopause 

2.  IR  radiance  approach.  In  this  technique  the  observed  IR  cloud  radiance 
is  used  to  deduce  a  cloud  top  temperature,  which  is  then  matched  to  a 
temperature  profile  specified  from  a  recent  analysis  or  short  term 
forecast. 

The  CDWs  simulated  by  Dey  ei  al.  essentially  used  an  error  free  temperature 
profile  for  height  assignment.  We  could  reassign  the  pressures  by  adding  an 
ad  hoc  error  to  the  CDW  pressures  directly,  or  by  assuming  an  ad  hoc 
temperature  error  and  reassigning  the  pressure  by  comparing  the  "observed " 
temperature  to  the  nature  run.  These  errors  would  have  to  be  fairly  large. 
These  approaches  would  indirectly  create  horizontally  correlated  errors. 

We  have  chosen  to  more  directly  induce  correlated  errors  in  the  CDW  data  by 
spatially  filtering  the  errors  in  the  simulated  data  We  divided  the  data  into 
high  and  low  winds  and  filtered  the  u  and  v  wind  components  for  each  group 
separately.  As  seen  in  figure  A.l  the  distribution  of  observing  pressures  has 
two  major  modes  and  we  choose  600  hPa  as  the  dividing  pressure  The  data 
are  for  12  UTC,  21  November  1979;  there  were  a  total  of  1935  CDWs.  The  filter 
we  use  is  a  weighted  average  followed  by  an  amplification  step.  The  averaging 
reduces  the  amplitudes  of  the  original  errors,  while  inducing  correlations  The 
amplification  step  multiplies  all  weighted  average  errors  by  a  constant  factor 
to  restore  the  original  amplitude  of  the  error  variance.  This  constant  is  equal 
to  the  variance  of  the  original  errors  divided  by  the  variance  of  the  filtered 
errors  The  weighted  average  operator  is  equivalent  to  the  first  pass  of  a 
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Barnes  filter  (Koch  et  al.,  1983)  The  weights,  w,  are  given  by  in  terms  of  the 
distances,  d,  between  the  point  to  be  filtered  and  the  nearby  points,  according 
to 


( 

w-exp 

V 


2\ 


(Al) 


Here  do  is  the  scale  of  the  filter  and  all  points  within  2xdo  are  included  in  the 
weighted  average.  After  filtering  a  random  error  field,  correlations  drop  to 
approximately  0.5  at  a  separation  distance  of  do.  We  chose  do  to  be  0.1  earth 
radii  or  approximately  637  km 


Fig.  A.l.  Distribution  of  CDW  pressures. 

The  results  of  the  filter  on  the  horizontal  correlations  of  the  CDW  errors  is 
shown  in  figures  A.2and  A.3  for  the  low  and  high  wind  groups.  For  reference 
the  correlation  of  the  nature  run  evaluated  at  the  same  data  locations  is 
presented  in  the  upper  row  of  plots.  Correlations  of  the  original  errors  are 
displayed  in  the  center  row  and  correlations  of  the  filtered  errors  are  displayed 
in  the  lower  row  of  plots.  As  indicated  the  results  for  the  u  and  v  wind 
components  are  displayed  in  the  left  and  right  columns  respectively.  Note 
that  the  errors  are  treated  by  the  filter  as  if  they  were  all  at  a  single  level 
(separately  for  high  and  low  winds).  On  the  other  hand,  the  nature  data  used 
for  the  correlations  are  at  a  variety  of  pressure  levels.  Tlie  synoptic  scale  is 
clearly  evident  in  the  v  correlations  from  the  nature  run.  This  is  not  obvious 
in  the  u  correlations  because  of  the  strong  zonal  component,  especially  at 
higher  levels. 
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Filtered  errors  Errors  Background 


Distance  (Earth  radii)  Distance  (Earth  radii) 

Figure  A.2:  Horizontal  correlations  for  low  level  (1000-600  hPa)  CDWs.  Shown  are 
correlations  for  the  background  field  evaluated  at  the  CDW  locations  (top),  original 
CDW  errors  (middle)  and  filtered  CDW  errors  (bottom)  for  u  (left)  and  v  (right) 
wind  components. 
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Filtered  errors 


0.08  0.16  0.24 

Distance  (Earth  radii) 


Figure  A. 3:  Horizontal  correlations  for  high 
A.2. 


level  (600-100  hPa)  CDWs.  Same  as  Fig. 


B.l  Methodology 


The  radiance  retrieval  is  a  step  in  CASCADE  following  the  OI,  but  preceding 
the  OI  jjostprocessor  (OIPTS).  As  implemented  in  our  OSSEs,  gridded  output 
from  the  OI  was  used  as  the  first  guess  for  the  retrieval,  and  updated  values 
were  used  at  grid  points  with  SSM/T  data.  The  horizontal  interpxslation  from 
analysis  grid  points  to  sounding  locations  (and  vice  versa)  was  not  explicitly 
simulated.  To  eliminate  difficulties  with  variable  surface  emissivities  which 
occur  over  land  or  ice,  retrievals  were  computed  only  at  grid  points  over 
water,  and  equatorward  of  65°.  Satellite  track  computations  were  used  to 
determine  the  data  swaths  for  each  6-hour  interval,  and  only  gridp)oint 
within  those  swaths  were  used  in  the  retrievals.  Out  of  the  total  of  6642 
analysis  grid  points  over  the  globe,  1782  (27%)  are  poleward  of  65°,  and  1402 
(21%)  of  the  remainder  are  over  land.  For  06  UTC  analysis  time  periods,  2269 
(34%)  of  the  remaining  points  were  not  covered  by  the  DMSP  satellite,  and 
1189  (18%)  were  used  in  the  retrievals.  Similar  distributions  of  points  apply  to 
the  other  analysis  time  periods. 

B.1.1  Retrieval  method 

The  formulation  of  the  simulated  errors  is  based  on  Eyre  (1988)  and  Eyre 
(1989)  analysis,  as  follows. 

Linear  optimal  retrieval,  like  many  other  retrievals,  may  be  written  as: 

^-^.  =  W(y«-yc(^o))/  (Bl) 

where  x  =  retrieval 
jc,  =  first  guess 
W  =  weight  matrix 
ym  =  measured  radiances 
yc  =  calculated  radiances 

Nonlinear  optimal  retrieval  involves  an  iteration 

=  ^,-1  +  W(y„  -  y,(x„_,))  +  B(j:„  -  Ar„.,)  (B2) 

where  B  is  a  matrix  controlling  how  much  and  in  what  way  the  retrieval  may 
deviate  from  the  first  guess. 

At  convergence  (i„  =  x„_■^  =  x),  so  that 
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(B3) 


i-x,  =  B-’W(y.-y,(i)) 

In  both  cases  we  have 

i-x,  =  B-'W(y, -y,(x.))  (B4) 

where  Xb  is  Xg  in  the  linear  case  and  x  in  the  nonlinear  case  and  B=I  for  the 
linear  case.  For  our  error  analysis  we  linearize  ym-  First 

ym=yc(xT)  +  £,+e,  (B5) 

where  Xj  =  true  state  vector 

=  forward  problem  errors  including  factors  left  out  of  x 
=  instrument  errors 


Then, 


yc(^r)=y.(j^B)  +  K(jr7--Xfc)  +  e,  (B6) 

where  K  =  derivative  of  forward  problem  with  respect  to  x,  evaluated  at  xj, 
£f  =  linearization  error 

Combining  the  last  3  equations: 

x-x^  =  x,-Xt+  B-’W(K(x^  -  ^»)  +  (B7) 

where  e^  =  +€,+ 

So  far  we  have  made  no  approximations;  modeling  will  require  several.  In 
the  linear  case,  x^  =  x„,  B  =  I  and 

x-j:r  =  (I-WKXx<,-Jr7-)  +  We„ 

For  the  linear  optimal  retrieval 

W  =  (C-’  +  K^’E-'K)''  K^E-’  =  CK’'(KCK^  +  E)"'  (B9) 

Here  E  is  the  covariance  matrix  of  the  measurement  and  forward  model 
errors,  C  is  the  covariance  matrix  of  the  background  errors,  and  K  is  the 
matrix  of  the  derivatives  of  y  with  respect  to  x,  that  is  the  derivative  of  the 
forward  problem. 

In  the  nonlinear  case,  x^  =  x  and 
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(BIO) 


(I  +  B-'WK)(i  -  JTr )  =  (x.  - 1, )  +  B-' W£. 
For  the  optimal  retrieval 


W  =  (k’'E-’K  +  C-’)’’k^E 

(BID 

B  =  (K^E-’K  +  C-’)'’c-’ 

(B12) 

Therefore 

(i-X7.)  =  B(x,-x^)  +  We„  (3^3) 

This  is  very  much  like  the  original  nonlinear  iteration  of  Eq.  (B2). 

To  apply  these  relationships,  in  general  we  need  W,  K  and  B  as  well  as  Xg  and 
xj-  and  a  way  to  model  the  errors.  For  optimal  retrievals  we  can  generate  W,  K 
and  B  from  K,  E  and  C  used  in  the  retrieval.  In  such  a  case  the  same  K  is  used 
in  the  retrieval  and  error  analysis.  This  is  permissible  so  long  as  it  is 
accounted  for  in  generating  that  part  of  the  error  due  to  the  linearization. 

The  definition  of  Xj  is  troublesome.  In  many  OSSE  setups  xj  does  not  contain 
enough  small  scales.  An  increment  should  be  added  to  the  nature  run  to 
account  for  scales  between  the  scales  observed  by  the  instrument  (the  volume 
of  the  atmosphere  it  averages  over  for  space  based  instruments)  and  the  scales 
represented  in  the  nature  run.  It  is  probably  safe  to  consider  these  increments 
independent  and  uncorrelated  from  each  other. 

Generating  the  errors  is  simple  if  they  can  be  modeled  as  normally  distributed 
random  variables  with  a  known  bias  and  channel-channel  covariance  matrix 
and  no  correlations  in  the  horizontal  direction.  These  statistics  might  even  be 
stratified  by  location,  surface  type,  etc.  For  example,  suppose  e  is  the  error,  b  its 
bias  and  E  its  covariance  matrix.  Then  let  A  be  the  matrix  of  eigenvectors  of  E 
and  D  the  diagonal  matrix  of  the  square  root  of  the  eigenvalues  of  E  (i.e. 
E=ADDA^).  The  modeled  error  is  then  given  by  b+ADu,  where  u  is  normally 
distributed  with  zero  mean  and  unit  variance. 

Of  course  the  errors  are  horizontally  correlated.  Even  if  the  instrument  error 
is  not,  the  forward  problem  and  linearization  errors  are.  If  these  horizontal 
correlations  are  small  compared  to  that  of  the  first  guess  error  iXg-Xf)  it  is  fair 
to  ignore  them.  To  assess  this  we  may  compare  estimates  of  both  correlations. 
The  magnitude  of  the  first  guess  error  correlations  are  fairly  well  known  from 
previous  studies  for  OI  schemes.  The  forward  problem  errors  might  be 
studied  in  simulation  using  a  more  accurate  radiative  transfer  model. 
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B.1.2  Pre-  and  postprocessii^ 

The  forward  problem  for  the  computation  of  the  radiances  uses  a  perfect 
input  profile  (bilinearly  interpolated  from  1°  by  1®  nature  data)  at  the  nature 
run  mandatory  pressure  levels,  with  errors  of  representativity  added  for  a  125 
km  averaging  length  scale.  The  standard  deviations  of  these  errors  are  0.7  K 
for  temfjerature,  and  1.3  m/s  for  wind  components,  and  0.1  for  relative 
humidity.  The  values  for  temperature  and  winds  are  based  on  Nastorm  and 
Gage.  The  value  for  humidity  is  ad  hoc.  Appropriate  instrument  errors  were 
added  to  the  radiances  computed  from  the  forward  problem. 

The  background  for  the  retrieval  is  obtained  from  the  analysis.  The  output 
from  the  OI  consists  of  analysis  values  of  temperature  T  and  specific  humidity 
q,  along  with  estimates  of  the  analysis  error  standard  deviations  of 
geopotential  height  z  and  relative  humidity  rh,  all  on  sigma  layers.  The  Eyre 
retrieval,  on  the  other  hand,  requires  background  values,  and  the  complete 
error  covariance  matrix  of  T  and  log  (q)  on  40  (15  for  moisture)  pressure  levels 
extending  to  0.1  hPa.  Thus,  a  certain  amount  of  pre-  and  postprocessing  was 
required. 

As  a  first  step,  the  analysis  values  were  used  to  compute  an  updated  surface 
pressure,  and  sigma  layer  pressures  were  derived  from  the  result.  The 
analysis  values  of  T  and  q  (converted  to  log(q))  were  inter /extrapolated  to  the 
mandatory  pressure  levels  (constant  values  were  used  for  extrapolation), 
along  with  the  estimated  analysis  errors  of  rh.  Interpolated  values  of 
moisture  were  checked  for  supersaturation. 

The  estimated  analysis  errors  of  height  (on  a  layers)  were  converted  to  those 
of  temperature  by  first  constructing  a  height  error  covariance  matrix,  using 
the  following  specified  error  correlation  structure: 

r«= - ^-p-  (B14) 

with  Vj„=1.098.  The  value  of  the  constant  was  obtained  by  fitting  this  curve  to 
the  RAOB  error  correlation  values  used  in  the  simulation  of  RAOB  data.  The 
covariances  of  height  are  then  related  to  those  of  temf)erature  through  (Eyre, 
1989): 

S^  =  HS,H^,  (B15) 

where  H  =  L'l,  L  being  the  operator  that  relates  T  to  z  (z  =  L.T),  and  S^,  S,  are 
the  error  covariance  matrices  of  heights  (z)  and  temperature  (T).  The  above 
relationship  follows  directly  from  the  definition  of  Sj  and  S,: 
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S,  =  <z  z^  >, 

Sj  =  <T  T'^  >  =  <{H  z)  {H  zV  >  =  H  <zz^  >H^, 


(B16) 


(BIT) 


where  <  >  denotes  the  expectation  operator. 

The  above  calculation  was  only  done  for  a  subset  of  the  a  layers  (the  bottom¬ 
most  and  top-most  layers  were  excluded),  because  otherwise  an  ill- 
conditioned  matrix  resulted  in  unreasonably  large  temperature  errors  (20  -  30 
K).  Using  only  the  central  o  layers  more  reasonable  values  for  the 
temperature  errors  (1  -  5  K)  were  obtained.  Temp)erature  errors  were  further 
constrained  explicitly  to  the  interval  0  to  4K.  Temperature  errors  for  the 
bottom  and  top  layers  were  set  equal  to  the  next  computed  layer.  We  removed 
the  barotropic  error  variance  from  the  height  errors  tefore  this  calculation, 
assuming  it  was  equal  to  the  height  error  at  the  bottom-most  o  layer.  The 
resulting  temperature  error  standard  deviations  were  then  interpolated  to  the 
mandatory  pressure  levels,  along  with  the  estimated  analysis  errors  of  rh,  but 
no  extrapolation  is  done  at  the  top.  Analysis  errors  of  rh  were  then  converted 
to  those  of  log(<7)  following  Eyre  (1989,  App.  B): 

lnq  =  \nq^  +  lnrh  (B18), 

which  implies 


d]nq  =  d\nq,+d\nrh  =  c^dT  (B19) 

with  c,  =  =  .08  K  ’,  and  c^=rh  =  .7,  so  that 

dT 

<ln^(y>=^-^^^-^  +  2^<r/iT>  +  c?<T^  >,  (B20) 

Cj  Cj 


and 

<rh^>  =  C2[<  Ini^^  >  -2  c,  <  In^  T  >  +  c?  <  >].  (B21) 

We  used  equation(B20),  with  the  additional  assumption  that  <rhT>  =  0,  to 
convert  <rh^  >  into  <lT\^q>,  and  (B21)  for  the  reverse  computation. 

All  mandatory  pressure  level  data  were  then  inter/extrapolated  to  the 
retrieval  pressure  levels  (using  the  same  matrices  as  in  the  simulation  of  the 
statistical  retrievals  of  TOYS  data).  Using  specified  error  correlation  structures 
(from  ECMWF  first  guess  error  statistics),  the  retrieval  level  error  covariance 
matrix  was  then  constructed,  and  the  Eyre  retrieval  computed. 
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The  output  of  the  Eyre  retrieval  was  then  interpolated  back  to  sigma  layers  as 
follows:  the  temperature  error  covariance  matrix,  and  the  error  standard 
deviation  of  rh  Ull  at  retrieval  pressure  levels)  were  extracted  from  the 
updated  error  covariance  matrix.  The  temperature  error  covariances  (on 
retrieval  p  levels)  were  then  converted  to  height  errors  (using  S,  =  L  •  Sj-  •  L), 
and  all  pressure  level  values  interpolated  back  to  sigma  layers.  Sigma  layer 
values  of  temperature  and  q  were  updated  by  the  interpolated  increment  of 
temperature  and  log(q),  estimated  analysis  errors  were  replaced  by  their 
interpolated  values  (if  they  were  smaller  than  the  input;  in  the  case  of  height, 
the  barotropic  component  was  added  back  in).  All  retrieved  and  interpolated 
values  were  checked  for  supersaturation. 

B.2  Retrieval  results 

Before  implementation  of  the  radiance  retrieval  in  the  OSSE,  some  off-line 
tests  were  conducted  using  the  analysis  from  the  CONTROL  OSSE  for  00  UTC 
24  Jan.  as  the  background  field.  These  tests  all  were  based  on  comparisons  of 
data  at  mandatory  pressure  levels  before  and  after  the  retrieval  steps. 

Examination  of  individual  profiles  showed  retrieved  values  closer  to  nature 
than  the  first  guess  at  most  levels,  while  some  levels  showed  increments  in 
the  wrong  direction,  and  some  in  the  right  direction,  but  overshooting.  It  was 
verified  that  all  the  interpolations  are  correct,  and  the  output  errors  of  z  and 
rh  are  properly  deaeased.  In  general,  the  actual  temperature  errors  of  the  first 
guess  are  smaller  (order  1-2  K  in  the  troposphere)  than  the  error  estimates 
(most  likely  an  artifact  of  the  height  to  temperature  conversion),  but  they 
nevertheless  show  the  correct  dependence  on  the  estimated  height  analysis 
errors,  and  the  output  height  errors  are  properly  decreased  due  to  the 
retrievals  (the  same  holds  for  the  rh  errors). 

Using  zero  instrument  and  representativity  errors,  retrieval  increments  of 
temperature  were  generally  small  in  the  troposphere  (rms  decrease  of  error 
compared  to  first  guess  of  order  .1  K  or  less).  At  the  tropopause  (150  hPa),  and 
at  300  and  250  hPa  there  was  a  slight  degradation  (less  than  .2  K).  Relative 
humidity  error  showed  a  uniform  reduction  due  to  the  retrievals. 

Adding  instrument  and  representativity  errors  degrades  the  retrievals 
slightly,  to  the  point  where  there  is  negligible  positive  impact  or  a  slight 
degradation  in  the  troposphere.  In  the  stratosphere  there  is  a  positive  impact 
(.08  K  to  1.27  K  between  100  hPa  and  50  hPa,  3.63  K  at  30hPa). 

Mandatory  level  error  statistics  were  also  collected  during  the  radiance 
retrieval  step  during  the  SSMTRAD  OSSE.  Fig.  B.l  shows  the  difference  in 
mean  and  rms  errors  of  the  analysis  grid  point  values  before  and  after  the 
radiance  retrieval  from  the  whole  seven  day  assimilation,  computed  only  for 
those  points  at  which  retrievals  were  performed.  In  general  terms,  these 
results  agreed  with  our  previous  off-line  calculations. 
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Errors  of  relative  humidity  show  a  uniform  improvement  between  850  hPa 
and  300  hPa,  with  reductions  up  to  20  %.  At  1000  hPa,  there  is  a  slight 
degradation.  The  rh  error  reductions  are  due  to  both  a  reduction  in  the  bias 
and  standard  deviation. 

Temperature  errors  show  almost  no  impact  in  the  troposphere,  except  that  a 
warm  bias  is  introduced  by  the  retrieval  at  250  hPa.  In  the  stratosphere,  there 
is  a  reduction  of  error  (on  the  order  of  0.5  K),  except  at  70  hPa  the  retrieval 
also  introduces  a  warm  bias,  which  leads  to  an  overall  degradation  at  that 
level.  At  the  topmost  mandatory  level  (10  hPa,  not  shown)  the  analysis  errors 
before  the  retrieval  are  quite  large  (16  K),  and  the  retrieval  impact  is  on  the 
order  of  10  K.  However,  since  the  topmost  o  layer  is  at  approximately  20  hPa, 
the  analysis  errors  are  due  largely  to  extrapolation  errors,  and  the  retrieval 
impact  is  largely  lost  in  the  interpolation  from  pressure  back  to  a  layers. 

The  small  impact  of  the  retrieval  on  temperature  errors  in  the  troposphere  is 
consistent  with  the  error  characteristics  of  the  instrument  and  the  small  first 
guess  errors.  However,  the  degradation  at  250  hPa  and  70  hPa  is  somewhat 
puzzling.  One  possible  explanation  is  that  the  error  covariance  matrix  used  in 
the  retrieval  was  incorrect,  a  distinct  possibility  since  it  was  derived  through  a 
number  of  steps  with  several  ad  hoc  assumptions  about  underlying  error 
characteristics.  Another  possibility  is  that  because  the  problem  is  so  severely 
underdetermined  (only  1 1  channels  were  used  in  the  retrieval  of  40 
temperature  and  15  moisture  values)  the  matrix  inversions  were  ill- 
conditioned. 
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The  moist  initialization  was  designed  to  explore  and  test  methods  of  using 
observed  precipitation  rates  in  the  CASCADE  data  assimilation  scheme.  In 
this  exploratory  study,  no  errors  were  added  to  the  nature  run  precipitation. 
The  impacts  thus  represent  an  upper  limit  for  the  techniques  employed  here. 
The  methods  were  implemented  as  an  initialization  step  and  a  forecast 
modification  of  the  first  guess  forecast.  The  forecast  modification  consisted  of 
a  rescaling  of  the  predicted  convective  heating  (and  precipitation)  by  the  ratio 
of  observed  to  predicted  convective  precipitation  for  all  grid  points  in  the 
tropics  with  non  zero  predicted  convective  rain.  The  factor  was  not  allowed  to 
exceed  the  maximum  value  of  10.  The  initialization  step  followed  all  other 
analysis /retrieval  steps  of  CASCADE.  In  the  MOIST  INIT  OSSE,  it  was 
inserted  after  the  OI  and  before  the  next  6-hour  forecast;  in  the  CASCADE 
OSSE,  it  followed  the  SSM/T  radiance  retrieval  step.  The  three  components  of 
the  initialization  step  consisted  of  a  diabatic  NMI,  a  separate  divergence 
adjustment  step  in  which  vertical  profiles  of  horizontal  divergence  were 
adjusted  based  on  observed  precipitation  rates,  and  a  temperature  and 
moisture  adjustment  step  designed  to  bring  predicted  convective 
precipitation  rates  in  closer  agreement  with  observations.  In  the  following, 
we  describe  the  three  components  and  their  combination  into  one 
initialization  procedure,  along  with  test  results  from  a  case  study  forecast. 
Results  of  the  moist  initialization  OSSEs  are  discussed  in  the  context  of  the 
other  OSSEs. 

Cl  Diabatic  NMI 

Implementation  of  the  diabatic  NMI  required  only  minor  changes  to  the 
global  spectral  model  (GSM).  Changes  were  made  to  include  physical 
processes  in  the  tendency  calculations,  and  to  convert  the  adjustments  in  the 
moist  physics  (Kuo  convection,  large-scale  precipitation,  and  dry  adiabatic 
adjustment)  into  tendencies.  In  the  tropics,  we  scaled  the  predicted  convective 
heating  by  the  ratio  of  observed  to  predicted  precipitation  (as  in  the  forecast 
modification  for  the  first  guess  forecast).  This  amounts  to  using  observed 
column  integrated  latent  heating  rates,  but  only  at  model  grid  points  with 
convection. 

The  same  2-iteration  Machenhauer  initialization  scheme  was  used  as  in  the 
adiabatic  NMI  in  the  other  OSSEs,  using  six  vertical  modes  and  a  48  hour 
frequency  cutoff. 

C2  Divergence  adjustment 
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C2.1  Determination  of  divergence  profile 

Julian  (1984)  has  shown  that  the  vertical  profile  of  divergence,  or, 
equivalently,  vertical  velocity,  can  be  related  to  cloud  top  equivalent  black 
body  temperature  (EBBT)  in  the  tropics,  a  quantity  closely  related  to 
convective  precipitation  rate.  Kasahara  et  al.  (1988)  devised  a  scheme  in 
which  a  parabolic  'ertical  profile  of  vertical  velocity  was  prescribed  for  points 
with  EBBTs  colder  than  a  threshold  value  (presumed  to  be  precipitating),  and 
the  maximum  value  was  linearly  related  to  the  difference  between  observed 
and  threshold  EBBT.  For  cloud-free  points,  vertical  velocity  was  determined 
as  a  residual  from  the  thermodynamic  equation. 

We  investigated  several  different  approaches  to  specifying  divergence  in  the 
tropics  from  the  observed  precipitation.  Motivated  by  the  approach  of 
Kasahara,  which  implies  a  linear  vertical  profile  of  divergence  for 
precipitating  points,  we  computed  linear  regression  statistics  of  divergence  (D) 
or  divergence  divided  by  convective  rain  rate  (D/R)  against  pressure.  Nature 
run  divergence  D  and  D/R  on  either  the  R40  transform  grid  (128x102)  or  T106 
grid  (320x160)  were  regressed  against  pressure  for  p>oints  with  R>e,  where  e  is 
a  specified  lower  limit.  The  convective  rain  rate  from  the  nature  data  was 
computed  as  a  6-hour  average  centered  around  the  analysis/initialization 
time,  and  interpolated  horizontally  as  necessary.  Results  from  a  test  of  this 
method  (using  data  from  12  UTC  January  20)  revealed  several  shortcomings, 
which  were  not  appreciably  different  for  the  R40  and  T106  grids.  Most 
importantly,  there  was  a  large  amount  of  scatter,  particularly  if  e  was  chosen 
small  enough  so  that  a  substantial  number  of  grid  points  were  included  in  the 
sample.  Correlation  coefficients  were  less  than  0.1  if  all  rainy  points  (53%  of 
all  points)  were  included,  and  less  than  0.2  if  only  31%  of  all  points  were 
included.  In  addition,  correlation  coefficients  were  larger  (although  only 
slightly)  for  D  than  for  D/R,  indicating  a  fundamental  problem  of  relating  the 
strength  of  the  ascent  to  the  intensity  of  precipitation  in  the  straightforward 
manner  attempted  here.  This  is  due  to  the  fact  that  R  varies  over  several 
orders  of  magnitude,  thus  making  the  quantity  D/R  a  highly  variable 
quantity. 

To  avoid  the  problems  encountered  in  the  linear  regression  techniques,  we 
computed  mean  profiles  of  vertical  velocity  and  divergence  for  several 
categories  of  grid  points  in  the  tropics,  each  category  being  defined  by  upper 
and  lower  limits  of  observed  convective  precipitation  rates.  This  was  done  for 
both  the  R40  and  T106  grids,  and  for  D  and  vertical  velocity  on  pressure  levels 
and  D  on  o  layers  (from  a  preprocessed  version  of  the  nature  run  data),  and 
using  both  convective  and  total  precipitation.  An  example  of  the  results  is 
shown  in  Figure  C.l,  which  shows  the  mean  profiles  of  D  (along  a  surfaces) 
for  5  categories  on  the  R40  transform  grid.  The  definition  of  the  categories  in 
terms  of  total  rain  rate  is  given  below.  Only  grid  points  equatorward  of  30° 
were  used  in  the  statistics. 
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Table  C.l;  Definition  of  categories  1  -  5.  Riow  “nd  Rupp  are  the  lower  and  upper  limits  of  rain 
rate  (note  that  10'^  mjs  corresponds  to  0.864  mmiday). 


Category 

Rlow 

(lO^m/s) 

Rupp 

(10-8m/s) 

Number  of 
grid  points 

Percentage  of 
grid  p>oints 

1 

- 

1 

17669 

68% 

2 

1 

5 

4026 

15% 

3 

5 

10 

1691 

6% 

4 

10 

50 

2496 

10% 

5 

50 

- 

230 

1% 

The  profiles  shown  in  Fig.  C.l,  which  were  based  on  3  days  of  nature  run  data 
(12  UTC  20  January  through  00  UTC  23  January)  show  a  clear  connection 
between  rain  rate  and  low  level  (upper  level)  convergence  (divergence). 

There  is  still  a  substantial  amount  of  scatter:  all  the  mean  profiles  except  for 
category  5  are  within  one  standard  deviation  of  one  another;  however,  the 
mean  profiles  from  the  aggregate  statistics  of  the  3  days  (6  time  periods) 
shown  in  Fig.  C.l  agree  well  with  those  computed  for  individual  time 
periods. 

We  tested  the  sensitivity  of  the  results  to  a  number  of  factors,  but  all  showed  a 
negligible  effect.  The  results  shown  in  Figure  C.l  are  virtually  identical  when 
convective  rather  than  total  rain  rate  is  used.  Profiles  for  divergence  along 
pressure  surfaces  (rather  than  o)  are  also  very  similar,  except  that  the 
differences  between  the  category  means  are  slightly  larger;  however,  standard 
deviations  at  low  levels  are  also  larger.  Profiles  of  vertical  velocity  are 
consistent  with  those  of  divergence,  showing  ascent  (descent)  for  rainy  (dry) 
categories.  Results  are  virtually  identical  for  the  R40  and  T106  grids. 
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Fig  C.l:  Profiles  of  divergence  along  a  surfaces  for  5  different  categories  of 
rainfall  rate  for  grid  points  in  the  tropics.  Category  2  is  shown  in  solid,  2  in 
dotted,  3  in  dashed,  4  in  dash-dotted,  and  5  in  dash-dot-dot-dot. 


C.22  Adjustment  procedure 

Based  on  the  tests  described  above,  the  following  divergence  adjustment 
procedure  was  implemented.  For  each  gridpoint  equatorward  of  30°,  the 
category  (1-5)  was  determined  from  the  observed  total  rain  rate.  The  analysis 
divergence  profile  was  then  replaced  by  the  divergence  profile  corresponding 
to  the  rainfall  category,  with  the  following  modification:  (1)  to  avoid 
unreasonably  sharp  gradients  at  the  boundary  of  the  adjusted  region,  analysis 
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and  specified  divergence  profiles  were  combined  as  a  weighted  average  at  the 
latitude  circles  at  the  boundary;  (2)  to  avoid  changing  the  vertical  average  of 
divergence,  the  vertically  averaged  difference  between  analysis  and  sp)ecified 
divergence  is  added  to  the  sf)ecified  divergence  profile.  Thus,  if  we  denote  the 
original  divergence  by  d  and  the  divergence  profile  of  category  i  by  c^  and 
their  respective  vertical  integrals  by  D  and  then  the  modified  divergence 
after  the  first  step  of  the  adjustment  (d^)  for  a  grid  point  with  category  i  is 
given  by: 


d^  =  (1-w)  d  +  w  (c>  +  D  -  CO  , 


(Cl) 


where  w  is  the  blending  weight  (w=l  for  all  latitudes  equatorward  of  30°,  and 
w=0.5  for  the  latitude  circle  poleward  of  30°).  The  blending  of  analysis  and 
specified  divergence  fields  at  the  edge  of  the  tropics  is  similar  to  the  procedure 
employed  by  Kasahara  et  al.  (1988).  To  avoid  changes  to  the  divergence  field 
outside  the  tropics,  a  technique  similar  to  one  employed  by  Kasahara  et  al. 
(1988)  was  implemented  here.  Rain-free  grid  points  (category  1  in  our  case) 
were  further  modified  in  a  second  step  by  subtracting,  at  each  level,  the 
horizontal  integral  of  the  changes  of  divergence  introduced  in  the  first  step, 
normalized  by  the  area  of  category  1: 


=  d'-5i,w 


w(d*  -d)dA 

J^^wdA 


(C2) 


Here  5,,  is  the  Kronecker  delta  (=1  for  points  with  category  1,  =0  otherwise), 
and  the  area  integral  extends  over  all  grid  points  with  non  zero  vo.  This 
ensures  that  the  horizontal  integral  of  divergence  in  the  tropics  is  unaffected 
by  the  divergence  adjustment. 

C3  Temperature  and  moisture  adjustment 

In  contrast  to  the  empirically  based  divergence  adjustment,  adjustment  of 
moisture  and  temperature  values  in  order  to  achieve  better  agreement  of 
predicted  and  observed  convective  precipitation  is  quite  model-dependent. 
With  the  present  version  of  the  model,  which  uses  the  Kuo  convection 
scheme,  we  could  draw  upon  the  experience  of  Donner  and  Rasch  (1989)  and 
Krishnamurti  et  al.  (1984).  In  Donner  and  Rasch,  the  initial  state  values  of 
moisture  are  adjusted  to  obtain  the  desired  precipitation,  and  temperature 
values  are  adjusted  to  obtain  the  desired  heating  rate.  The  adjustments  are 
performed  as  part  of  a  variational  scheme  under  the  additional  constraint  of 
minimizing  the  column  integrated  squared  adjustments.  Because  derivatives 
of  precipitation  and  heating  with  respect  to  the  input  profiles  of  temperature 
and  moisture  are  discontinuous  across  the  point  of  conditional  instability, 
low-level  values  of  temperature  and  moisture  have  to  be  adjusted  separately 
to  ensure  the  presence  of  convection,  and  these  values  are  not  included  in  the 
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variational  adjustment  process.  However,  even  then  there  will  be  some 
points  for  which  the  criteria  for  convection  will  not  be  satisfied,  particularly  if 
there  is  no  or  insufficient  convergence  of  water  vapor.  It  was  found  in 
Kasahara  et  al.  (1992)  that  the  separate  divergence  adjustment  step  was 
necessary  to  improve  the  results  of  the  temperature  and  moisture 
adjustment. 

Krishnamurti  et  al.  (1984)  use  a  somewhat  different  approach.  In  their  work 
they  use  a  "reverse  Kuo",  in  which  the  preliminary  predicted  (not  the  initial 
values)  of  moisture  are  adjusted  by  a  factor  constant  with  height  (1+e).  The 
adjustment  step  is  quite  straightforward  in  this  case,  since  the  amount  of 
moisture  convergence  can  be  directly  prescribed  by  the  proper  choice  of  e. 
However,  their  approach  requires  a  forward-backward  integration  during  an 
extended  pre-integration  period  that  is  as  long  as  2  days. 

We  used  a  combination  of  the  two  approaches  in  our  scheme.  As  in  Donner 
and  Rasch,  the  initial  state,  rather  than  the  preliminary  predicted  values  of 
temperature  and  moisture  were  considered  for  adjustment.  To  include  as 
many  points  as  possible  in  the  adjustment  procedure,  restrictions  on 
convection  (minimum  values  for  moisture  convergence)  were  removed,  and 
surface  values  of  moisture  and  temperature  were  increased  at  grid  points 
which  were  stable  but  had  non  zero  observed  precipitation.  Both  the  lowest 
level  temperature  and  specific  humidity  were  increased  by  a  predetermined 
increment  until  a  convective  cloud  deep  enough  for  the  production  of 
precipitation  was  diagnosed,  or  the  maximum  number  of  iterations  was 
exceeded. 


We  also  attempted  a  simplified  version  of  the  moisture  adjustment,  in  which 
the  initial  moisture  values  above  the  surface  were  modified  by  a  factor  (1+e), 
similar  to  Krishnamurti's  approach.  A  closed-form  solution  for  this  factor 
was  derived  from  an  approximation  of  the  modified  convective  precipitation 
resulting  from  this  adjustment.  This  derivation  is  reproduced  in  the 
following: 

In  the  version  of  the  Kuo  scheme  implemented  in  our  GSM,  the  convective 
precipitation  rate,  R,  is  given  by  (neglecting  the  reevaporation  of  falling 
precipitation) 


R  =  ^.w'  (l-b)  , 
g 


(C3) 


where  P,  is  the  surface  pressure,  g  the  acceleration  of  gravity,  w  the  moisture 
convergence,  and  (1-b)  the  heating/moistening  partitioning  parameter.  The 
moisture  convergence  is  given  by  the  vertical  integral 
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dt  ^ 


k-K 


vdtA 


A<Tk  . 


(C4) 


where  a,,  k,  denote  the  cloud  top,  and 


dq"^ 

^dt 


is  given  by 


yk 


^  K*'-qk')/(2At)  fort>2 
.‘^Wk  l(qL*'-ql)/A t  fort  =  1, 


since  a  forward  step  is  used  for  the  first  time  step,  and  leap-frog  steps 
thereafter.  Finally,  the  parameter  b  in  this  scheme  is  computed  from 


b  = 


TCp  dd 

3 

jO) 

d  L  dp 

dp 

dp 


2fca^dp 

J  dp 

fa 

with  (a  approximated  as  o  =  ir  =  -  I  V  Vd<T,  and  the  vertical  integrals 
approximated  by 


In 


b  = 


^  ®k)'*'(qk*i  9k) 


I 

en 


(C5) 


(C6) 


2X^k(qk*i  qk)  (C7) 

The  sums  are  evaluated  from  cloud  bottom  to  cloud  top. 

To  assess  the  effects  of  modifying  the  initial  moisture  on  the  rain  rate,  we 
need  to  consider  the  effects  on  w  and  1-b  separately.  The  moisture 
convergence  rate  may  be  written  as 

iy=  r^d(T=  f-V(Vq)da  -  [dq]^  +  r  T (CT)d(T+ fq(C  +  ^^^^) do,  (C8) 

Ja,  Ja,  *  "•  Ja,  J  * 

where  the  conservation  equation  for  moisture  was  used  to  substitute  for 

a  I 

T  (<j)  is  the  turbulent  tendency,  and  C  =  V •  Vlnp,.  If  we  now  consider  a 
modified  moisture. 


q“'=(l  +  e)q 

the  resulting  modified  moisture  convergence  is  given  by 


(C9) 
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(CIO) 


«;"'=(!  +  e)£  q(-V  •  V)  da  -  £  V  V((l + e)q) 

-Kl  +  e)[dql^  +  flT(a"’)  -  T(a)l  da  +  (1  +  £)£  q[C  +  ^-Inp.lda. 

If  we  express  the  turbulent  tendenqr  as  a  flux  divergence,  the  vertical  integral 
reduces  to  the  difference  between  the  turbulent  flux  at  cloud  top  and  the 
evaporation  at  the  surface.  We  may  safely  neglect  the  flux  at  the  top,  and, 
since  the  lowest  level  moisture  will  not  modified  according  to  (C9),  the 

turbulent  contribution  to  w"  -  w  can  be  neglected.  The  vertical  advection 
term  can  also  be  neglected,  since  d  =  0  at  a  =  1 ,  and  dq  is  small  at  a  =  a, . 
Further  neglecting  V(e  q)  leaves  the  approximate  relationship 

6„=w’"-w  =  ea^ 

with  a^  =  f  q(-V- V)  da+ f  q[C  +  ^lnpJ  da,  (Cll) 

which  can  be  evaluated  inside  the  Kuo  convection  subroutine. 

The  moistening  parameter  is  affected  indirectly  by  the  modification  of  the 
initial  moisture,  because  its  computation  is  based  on  the  (t+1)  time  level 
quantities.  If  we  neglect  effects  of  the  modified  q  at  time  t  on  the  values  of 

and  approximate 

q™'*' =(l  +  £)q'"‘*',  (C12) 

which  is  equivalent  to  assuming  that 


=  (l  +  e) 


St 


(C13) 


the  integral  involving  q  is  multiplied  by  the  factor  (1  +  e),  which  results  in  the 
approximate  formula 

=  (1  -  b)-"  -  (1  -  b)  =  [|  -  (1  -  b)]^.  (C14) 

2  l  +  £ 

Combining  the  changes  to  w  and  1-b  results  in  the  modified  rain  rate 

5  R  =  R"  -  R  =  Iw6, ,  +  (1  -  b)<5,  +  5,-b5w]—  (CIS) 

g 
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Setting  5  R  =  R°*"  -  R  and  denoting  Sr  =  -^  [R°**  -  R),  results  in  the  following 
equation  for  e  after  some  algebraic  manipulations: 

e^-  Ci  +  e  cj+Cj  =0. 

with  c,  =^a„,C2=(l-b)a„+[^-{l-b))  w-5„C3  =  -5r.  (C16) 


C3  Tests  of  the  initialization  procedure 

All  components  of  the  initialization  procedure  were  tested  on  the  last 
analysis  from  the  spinup  experiment  (00  UTC  23  January).  A  number  of 
different  configurations  were  tested,  in  order  to  evaluate  the  different 
initialization  components  in  isolation  and  in  combination.  The  following 
table  provides  a  summary  of  the  combination  of  initialization  steps.  The 
primary  tool  for  assessment  was  a  comparison  of  the  adjusted  initial  state, 
and  of  the  6-hour  forecast  generated  from  it,  with  the  corresponding  nature 
data. 

Table  C.2:  Configuration  of  initialization  experiments.  Convective  heating 
and  precipitation  scaled  by  ratio  of  observed  to  predicted  precipitation 
during  NMI  and  forecast  in  all  experiments  other  than  CONTROL 
Temperature  and  moisture  adjustment  refers  to  the  adjustment  of 
lowest  layer  values  only. 


Name 

Configuration 

CONTROL 

adiabatic  NMI 

Test  A 

diabatic  NMI 

Test  B 

Test  A  +  divergence  adjustment 

Teste 

Test  B  +  temperature  and  moisture  adjustment 

Test  D 

Test  A  +  temperature  and  moisture  adjustment 

Test  E 

Test  C  +  diabatic  NMI 

Test  F  (MOIST  INIT) 

Test  E  +  temperature  and  moisture  adjustment 

67 


The  diabatic  NMI  by  itself  (Test  A)  lead  to  very  small  differences  when 
compared  to  the  adiabatic  NMI  (CONTROL):  temperature  differences 
(vertically  averaged  rms  differences)  on  the  order  of  .09  K  and  wind 
differences  of  .07  m/s  in  the  initial  state,  and  rms  differences  from  nature 
agreed  to  within  .003  K  and  .001  m/s.  After  a  6-hr  forecast,  differences 
compared  to  CONTROL  were  still  small  (.1  K,  .3  m/s).  This  result  is  consistent 
with  the  findings  of  Kasahara  et  al.  (1992)  and  Knowlton  et  al.  (1989),  among 
others.  Interestingly,  forecast  errors  after  6  hours  were  worsened  for 
temperature  and  winds  by  the  diabatic  NMI  when  compared  to  control, 
whereas  the  moisture  errors  were  improved;  rms  errors  differed  by  only  .003 
K,  .02  m/s,  and  .01  g/kg,  however. 

Adding  the  divergence  adjustment  after  the  diabatic  NMI  (Test  B)  lead  to  the 
most  noticeable  changes.  In  the  initial  state,  the  vertically  averaged  rms  error 
of  the  divergent  wind  was  decreased  by  nearly  1  m/ s  (from  3.861  m/s  to  2.864 
m/s).  The  improvement  due  to  the  divergence  adjustment  is  also  reflected  in 
the  number  of  grid  points  with  both  observed  and  predicted  precipitation: 
during  the  first  time  step  of  the  forecast,  11.5  %  of  all  grid  points  had  both 
observed  and  predicted  rain  (hit),  40.9  %  had  observed  but  no  predicted  rain 
(miss),  and  3.7  %  had  no  observed  but  predicted  rain  (false  alarm)  after  the 
divergence  adjustment.  The  corresponding  numbers  for  the  diabatic  NMI 
alone  (Test  A)  were  10.4  %  hits,  42.0  %  misses,  and  4.5  %  false  alarms.  After  a 
6-hour  forecast,  rms  errors  of  temperature,  winds,  and  moisture  were  all 
decreased  compared  to  CONTROL  or  the  diabatic  NMI  alone  (Test  A),  but  the 
improvement  in  the  divergent  wind  errors  was  reduced  to  .4  m/s.  The 
improvement  in  the  hit/ miss  statistics  of  observed  and  predicted  rain  did  not 
persist  past  the  first  hour  into  the  forecast. 

The  output  from  the  divergence  adjustment  (Test  B)  was  subjected  to  the  first 
part  of  the  temperature  and  moisture  adjustment  (in  which  only  the  bottom 
layer  temperature  and  moisture  were  modified  as  necessary  for  conditional 
instability)  in  the  next  experiment  (Test  C).  In  the  initial  state,  this  improved 
the  lowest  level  specific  humidity  error  by  .13  g/kg  (from  2.76  to  2.63  g/kg); 
more  importantly,  the  hit/miss  statistics  at  the  first  time  step  of  the  forecast 
were  further  improved  to  13.2  %  hits,  39.2  %  misses,  and  3.8  %  false  alarms. 

As  in  the  case  of  the  divergence  adjustments,  however,  the  improvement  in 
the  hit/miss  statistics  did  not  persist  past  the  first  hour  of  the  forecast.  At  the 
end  of  the  6-hour  forecast,  rms  errors  of  temperature,  winds,  and  specific 
humidity  all  showed  improvements.  Compared  to  CONTROL,  errors  were 
decreased  by  .02  K  (temperature),  .4  m/s  (divergent  wind),  and  .02  g/kg 
(specific  humidity);  improvements  compared  to  the  divergence  adjustment 
forecast  (Test  B)  were  roughly  an  order  of  magnitude  smaller. 

The  temperature  and  moisture  adjustment  of  the  lowest  layer  was  also 
applied  to  the  output  of  the  diabatic  NMI  in  Test  D,  to  test  the  importance  of 
the  divergence  adjustment  as  a  first  step.  As  might  be  expected  from  the 
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comparatively  large  effect  of  the  divergence  adjustment  alone,  results  from 
this  experiment  were  noticeably  worse.  The  initial  hit/miss  statistics  were  11.6 
%  hits,  40.8  %  misses,  and  4.7  %  false  alarms,  while  forecast  errors  were  larger 
than  in  Test  C  for  all  flelds,  and  larger  than  even  CONTROL  and  Test  A  for 
temperature  and  winds.  It  is  interesting  to  note  that  the  temperature  and 
moisture  adjustment  actually  lead  to  a  slight  degradation  of  the  temperature 
and  wind  forecast,  whereas  in  Test  C  it  slightly  improved  those  errors. 

The  second  component  of  the  moisture  adjustment,  which  was  applied  to  the 
layers  above  the  surface,  did  not  achieve  the  desired  result.  It  increased  the 
initial  error  in  specific  humidity,  worsened  the  hit/miss  statistics  of  the  first 
time  step,  and  lead  to  increased  forecast  errors  of  specific  humidity  (even 
when  compared  to  control).  There  are  several  possible  reasons  for  this:  one  or 
more  of  the  approximations  made  in  the  derivation  of  the  equation  for  e 
could  be  invalid,  and  there  could  be  an  error  in  the  derivation  or  the  coding 
of  the  scheme.  Because  of  the  limited  level  of  effort  available  for  this  task, 
these  questions  could  not  be  pursued  further.  The  temperature  and  moisture 
adjustment  was  thus  restricted  to  the  adjustment  of  the  lowest  level 
temperature  and  moisture,  which  had  a  small,  but  positive  impact. 

Kasahara  et  al.  (1992)  experimented  with  several  different  combinations  of 
their  initialization  steps,  and  found  the  following  arrangement  to  give  the 
best  results  (in  terms  of  the  precipitation  spinup  behavior): 

step  1:  diabatic  NMI 

step  2:  divergence  adjustment 

step  3:  temperature  and  moisture  adjustment 

step  4:  diabatic  NMI 

step  5:  temp)erature  and  moisture  adjustment 

We  tested  adding  steps  4  (Test  E)  and  5  (Test  F)  to  the  output  from  our 
temperature  and  moisture  adjustment  and  also  found  the  best  results  when 
all  5  steps  were  used,  although  the  additional  improvements  were  quite 
small:  compared  to  Test  C,  the  initial  errors  were  reduced  by  .003  K 
(temperature),  .03  m/s  (divergent  winds),  and  .002  g/kg  (specific  humidity).  At 
the  end  of  the  6-hour  forecast,  rms  errors  were  decreased  by  .04  m/s  for  the 
divergent  wind,  while  other  quantities  showed  inconsistent  further 
improvements.  Even  though  the  additional  improvements  were  minor,  the 
moisture  initialization  used  in  the  MOIST  INIT  OSSE  consisted  of  all  5  steps 
listed  above  (Test  F  in  Table  C.2). 

The  characteristics  of  the  precipitation  spinup  are  depicted  in  Fig.  C.2  for  the 
various  experiments.  Shown  are  the  global  rate  of  convective  precipitation, 
stratiform  precipitation,  and  evaporation  for  each  time  step  during  the  6-hour 
forecast.  It  should  be  noted  that  the  convective  precipitation  rate  is  scaled  by 
the  ratio  of  observed  total  precipitation  to  predicted  convective  precipitation 
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in  the  tropics  in  all  experiments  other  than  CONTROL.  It  is  interesting  to 
note  that  in  our  model,  even  in  the  CONTROL  experiment,  there  is  very  little 
of  the  spinup  seen  in  other  models  (e.g.,  Kasahara  et  al.  1992).  Although  both 
convective  and  stratiform  precipitation  undergo  some  transient  behavior 
during  the  forecast,  initial  and  final  values  after  6  hours  differ  by  less  than 
30%  in  CONTROL.  The  initial  convective  precipitation  is  increased  by  the 
moist  initialization,  but  it  eventually  equilibrates  at  a  slightly  lower  level 
than  in  the  CONTROL  run.  Comparison  of  the  curves  for  Test  C  and  Test  F 
shows  that  the  additional  initialization  steps  reduced  the  amount  of 
overshooting  during  the  first  hour. 

The  stratiform  precipitation  shows  little  difference  between  the  exp)eriments 
during  the  first  hour,  but  a  clear  separation  of  the  exp>eriments  without 
divergence  adjustment  (CONTROL,  Test  A,  and  Test  D)  from  those  with 
divergence  adjustment  (Test  B,  C,  E,  and  F)  in  the  remainder  of  the  forecast 
pjeriod.  Final  values  of  the  divergence  adjustment  group  are  slightly  lower 
than  CONTROL.  The  evaporation  rate  at  the  beginning  of  the  forecast  is 
slightly  lowered  by  both  the  divergence  adjustment  and  the  temperature  and 
moisture  adjustment,  but  all  experiments  converge  to  nearly  the  same  value 
at  the  end  of  6  hours.  For  evaporation,  the  adjusted  exporiments  exhibit  more 
of  a  spinup  than  CONTROL.  This  indicates  one  of  the  weaknesses  in  our 
approach:  changes  in  the  divergence  and  low  level  temperature  and  moisture 
designed  to  yield  better  agreement  between  observed  and  predicted 
convection  also  affect  evaporation,  and,  indirectly,  the  moisture  convergence 
available  for  convection. 


/ 
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SDinup  during  6 -hour  lorecast 
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A  (dotted),  Test  B  (short  dashes),  Test  C 
(dash-dotted),  Test  D  (dash-dot-dot-dot).  Test  E 
(long  dashes),  and  Test  F  (solid). 


